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INTRODUCTION 


Three-dimensional  transonic  flows  through  a compressor  blade  row 
represent  a very  complex  problem.  Some  of  the  complications  that  make  the 
problem  very  difficult  to  solve  are  the  presence  of  adjacent  zones  of  subsonic 
and  supersonic  flow,  and  the  communication  between  these  zones  that  takes  place 
through  radial  pressure  gradients. 

As  a first  step  in  applying  recent  computational  techniques  to  this 

problem,  a study  of  the  nonlinear  small-disturbance  approximation  has  been 
fl-4) 

made  . This  report  contains  a description  of  the  computer  program  used  to 
obtain  the  numerical  results  presented  in  these  papers.  For  convenience,  the 
basic  equations  used  are  summarized  in  Section  1.  The  program  itself  is  des- 
cribed in  Section  2,  and  a sample  case  is  discussed  in  Section  3.  A listing 
of  the  program  and  a guide  to  preparing  the  input  are  given  in  the  Appendices. 

The  finite-difference  equations  used  are  the  set  on  which  the  results 
of  Reference  4 were  based.  They  supersede,  in  a few  small  details,  the  earlier 
sets  used  in  References  1-3. 

It  should  be  stressed  that  the  program  cannot  handle  cases  where  the 
inlet  relative  Mach  number  is  supersonic  at  the  tip,  because  the  farfield 
boundary  condition  used  is  not  a radiation  condition.  Thus,  Mach  waves,  origi- 
nating at  the  blade  row,  reflect  from  the  grid  boundaries  instead  of  escaping 
as  they  should.  As  a result,  the  solution  will  contain  a set  of  waves  which 
pass  back  and  forth  between  the  upstream  and  downstream  edges  of  the  grid. 

It  should  also  be  pointed  out  that  some  of  the  program  options  have  not 
been  used  extensively.  In  particular,  only  a few  runs  have  been  made  in  the 
design  mode,  where  the  loading  and  thickness  distributions  are  prescribed.  For 
this  reason,  there  are  unknown  limits  on  the  usable  ranges  for  some  of  the 
parameters  such  as  step  size  and  relaxation  factors. 
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Section  1 

SUMMARY  OF  BASIC  EQUATIONS 


The  coordinate  system  and  rotor  geometry  used  are  shown  in  Figure  1 , 
where  the  dimensionless  variables  are  defined  as 


(X>X  cJr 

^ = > p = 


U> 


> n:  = B-i  <j>  = ^ 


The  velocity  components  seen  by  a blade-fixed  observer  are 


K = -f-  U. 


, \AJf.  = V » 


\aIq  = o)r  + uT 


These  dimensional  perturbation  velocities  are  related  to  the  velocity 
potential  by 
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The  velocity  potential  satisfies  the  equation 

{j  - 4>Ai. 

+ i 1 


(1) 


(2) 


(3) 


(4) 


Subscripts  denote  derivatives  in  the  z,  p , coordinate  system;  thus  the 

symbol  (p^  stands  for  ^ • There  are  3 blades  in  the  row,  and  they 

are  taken  to  lie  in  the  helical  surfaces  defined  by  and  cJr  . The  axial 

projection  of  their  chord  is  a constant,  Cg^  . Thus,  they  are  located  at 


to  C-  2 'i  TT 

0 ^ Z.^  ; ^ ^ J— 


IXa 
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^ - 0 , 1 , X,. . B - 1 (5) 
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If  the  program  is  run  in  the  off-design  mode,  the  blade  shapes  are 
given  by  (see  Figure  2) 

ila  (S,  r)  = 4,(5. r)  t 1 i(s,r)  - oc(r)  ■ 5 (6) 

where  and  ot  are  the  camber,  thickness,  and  angle  of  incidence.  The 

blade-surface  boundary  condition  is 

K 


dn 

oLs 


(7) 


When  the  program  is  run  in  the  design  mode,  the  prescribed  quantities 
are  the  loading  and  thickness  distributions: 


AC^iS,r) 


z 


-Pu-  Pu. 

foo 


- Z (}  f P^) 


dn^.  oLf^^ 

ds  ds 


(8) 

(9) 


The  blade-surface  boundary  conditions  are  applied  in  the  helical  surfaces 
C = 0 and  ^ = zn/8  . 


Far  upstream  of  the  blades,  the  perturbation  potential  is  set  equal 
to  zero;  far  downstream,  it  is  expressed  as  <P  = C£  , where 
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n (i-M, 


-B 


(10) 
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Section  2 

PROGRAM  DESCRIPTION 


A guide  to  the  preparation  of  the  input  is  given  in  Appendix  A,  while 
the  listing  of  the  program  itself  and  a dictionary  of  the  FORTRAN  variables 
are  given  in  appendices  B and  C. 

The  program  begins  by  reading  and  printing  input  values.  The  finite- 
difference  grid  is  then  calculated,  in  subroutine  GRID.  The  region  in  which 
the  solution  is  to  be  found  is  divided  into  a grid,  with  the  indices  L , K 
and  N used  to  number  points  in  the  ^ , 2 and  p directions,  respectively. 

Equal  spacing  is  used  in  the  and  p-  directions: 

i;a)  = (L-  t)  , L = LMK  ■,  ^ /(LMX  - i)  (H) 

p(N)  = p^+(N-l)Ap  , A/  = Ap  = 


The  spacing  in  the  Z-  direction  is  nonuniform;  the  2-  coordinate  is  taken  as 


z = zcr; 


d ^ r J_  r = ^ 

dz  ^ dt  ’ ^ ' cLZ 


(13) 


The  variable  t is  then  allowed  to  vary  from  -1  to  +1,  as  H varied  from-®  to 

+ 0O  , 


r = -1  + KAZ 


K = t,  Z , i KMX  ; AT 


Z 

KMX  i-  1 


(14) 


The  particular  dependence  ZCT)  used  here  is 


where 
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The  locations  of  the  upstream  and  downstream  edges  of  the  grid  are  set  by  the 
input  values  for  FXB  and  FXI. 

Next,  the  blade-surface  boundary  conditions  are  calculated,  in  sub- 
routine BVAL.  These  conditions  depend  on  IBC:  for  IBC  = 1 (the  off-design 
case) , the  blade-shape  parameters  in  the  array  BV  are  used  to  calculate  the 
surface  slopes  dNDS  on  the  suction  and  pressure  sides,  at  each  of  the  axial 
and  radial!  grid  points  on  the  blades. 


For  IBC  = 2 (the  design  case)  the  assigned  loading  and  thickness 
distributions  are  used  to  generate  the  following  quantities,  which  are  used  in 
applying  blade-surface  boundary  conditions  at  L = 1 and  L = LMX,  respectively: 


DNDS (KB,N, 1) 

DNDS (KB,  N, 2) 


1 + 


ACp 

2 


-h  pi'(S) 


(17) 


(18) 


The  program  version  listed  here  uses  a parabolic-arc  blade  shape,  for 

IBC  = 1: 


- Ic"nr  (19) 

(20, 

The  maximum  thickness  and  camber,  which  occur  at  S/c  = 1/2  for  this  shape,  are 
allowed  to  vary  in  the  radial  direction  according  to: 

= C^\B\/(n  + BV(Z)  ^ + BV(3)  ^21) 
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VL. 


c Ibvc^)  + B \Ji5)  -^  + e /(<i) 

^ P p.  J 


(22) 


Note  that 


CiD 


= 1 / /r  + p^' 


(23) 


In  addition,  the  angle  of  incidence  varies  radially  as: 

a = B[/C7)  -!■  BV(S)  ^ + BVC9) 

P Pr 


(24) 


The  angle  of  incidence  is  measured  with  respect  to  the  helical  direction,  as 
shown  below: 


For  IBC  = 2 (the  design  case)  the  program  version  listed  here  uses 
a parabolic-arc  thickness  distribution: 

t\s)  = 4 -z  4-] 

c(r)  L C J 

_ ^ 1 

Co-  VTTp^ 


Z (K) 


(25) 
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The  maximum  thickness  is  taken  as  a constant; 


BV(3) 


(26) 


The  loading  distribution  had  the  form 


where  the  value  of  ACp^Cf)  varies  linearly  between  prescribed  values  at  the 
hub  and  tip,  and  a is  constant.  These  parameters  are  read  in  by: 

BV(Z)  = a , BVo)  = , Bm)  = (27) 

As  noted  in  Reference  3,  this  distribution  is  inconsistent  with  the  condition 
V = 0 at  the  hub  and  tip;  the  loading  A ()>(.!")  should  have  zero  radial  gradient 
at  the  hub  and  tip,  if  v is  to  equal  zero  there. 

Other  blade  shapes  and  loading  distributions  can  be  used,  by  making 
appropriate  changes  in  this  subroutine,  and  in  the  array  BV.  The  sole  function 
of  BVAL  is  to  return  the  arrays  DNDS  (KB,  N,  J),  J = 1,  2 as  defined  in 
Equations  (7),  (17)  and  (18)  above.  Use  of  other  blade-shape  or  loading  distri- 
butions does  not  require  changes  elsewhere  in  the  program. 
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Following  this  calculation  of  blade-surface  boundary  conditions,  the 
main  iteration  loop  begins.  The  outer  loop  (DO  300  N=...')  is  in  the  radial 
direction.  Within  this  loop,  the  line  on  which  the  solution  is  being  updated 
is  swept  downstream  in  the  loop  DO  5 K=..,  . The  innermost  loop  (DO  7 L=...') 
carries  out  the  relaxation  by  the  standard  recursive  algorithm,  described  below. 


The  finite-difference  equation  solved  in  the  innermost  loop  is 


m:.,  {c-i  (^t  - X - Vi  (t-.-  >:-j} 


-z 


(1+P  ) At  V i-,  A/  t-i  1 

I + <^/<  + ^ <t>^„  + </>^  ~ 


ZA  ^ 


O+pY/AT^ 

777  i 


-I 


'*P 


i%)  - 


(28) 


where 


'^K  " ’ - O +/>'')  - (il+1)  f -Aj<±2 ^ !<_-■> 

w K 


zat 


M = 

1 ) 


■For  V ^ O 

K 


(29) 
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This  equation  is  rewritten  as 


L HTL*f  L ^-u  Nl 


* 
l\ 

'X  ' K 


fit  t - j 


x 'Tk 


- D , /.  = 2,  3,  ... , ^/wx  -/ 


(30) 


The  solution  is  found  by 


--  NlL*f 


4>  = E£(L)  + FF(L) 

K 


(31) 


where  EE(L)  and  FF(L)  obey  the  recursion  relations 


££(1)  = / [s^  + ££U-0] 

= [ 0^  - FF(L-f)]  / [8^  i-  ££(l-r)]  (32) 


Subroutine  BVALO  uses  the  blade-surface  or  periodicity  conditions  to  set  EE(1) 
and  FF(1).  Then  the  recursion  formulas  above  are  used  to  find  EE(L)  and  FF(L) 
for  L = 2,  3,  LMX  - 1.  Subroutine  BVALI  then  sets  the  value  of  the 

potential  at  LMX,  by  again  using  the  blade-surface  or  periodicity  conditions. 
Equation  (31)  is  then  used  to  find  the  values  of  the  potential,  from  LMX  - 1 
down  to  L = 1.  These  tentative  values  of  (p  are  stored  in  the  array  SV;  the 
updated  values  are  found  by 


sya)  (t-u»^<t>‘- 

K 


(33) 


The  specific  formulas  used  in  subroutine  BVALO  are: 
(a)  Periodicity  condition,  for  Z < o 


EE  (0*0  , FF  (t)  = 


(34) 
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(b)  Periodicity  condition,  for  Z > 

EE  a)  ^ o ^ FF  (1)  = + A<^(N) 

K 

(c)  Blade-surface  condition  for  IBC  = 1:  the  derivative  (j>g^  is 
approximated  by 


(P  ) = I - ^ g 

a2^  I }Z  Al^ 


r - K - 7 (Pm  P cLn\ 

1 1ZA2;  2 d.s>(^, 


_ + 3C^);  - 4,  <p'^_,  + (p'^ 


zil+p^) 


toA^i; 


This  expression  is  then  used  to  write  the  potential  equation  at  L = 1 as 

{<  * /}  (' -A*;)  {f..i  ”C, -.f,  [ ^ >;] , 

* K<:,  {f«.x  (zx  - X - x.jj 

^FX{f,^r(X..-X)-  f.-zCX  - X.,)'j- 


4 (^J  T -T  ~i., 


P C ^ 

Zii+P^)  At 


/Xm  + ^ I 


^ \ \ \ 'V*>,  1 , / , AyO  > A/-7,  7 1 ' 1 

7;  I (’*  T-/* 
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This  is  then  written  as 


^ '^4>^  •*■  CE  o (38) 

which  leads  to 

E E U)  = - BE  / flE  , FF  O)  = - CE / RE 

(d)  Blade- surface  condition  for  IBC  = 2: 

EEa)  - , FFa)  = + f ^ d z 

The  latter  integral  is  calculated  in  BVAL  before  the  iterations  begin,  and  is 
stored  in  the  array  DNDS  (K-KLEP,  N,  1),  where  KLEP  is  the  K-value  of  the 
station  just  upstream  of  the  leading  edge. 


(39) 


(40) 


The  specific  formulas  used  in  subroutine  BVALI  are: 

(a)  Periodicity  condition,  forZ<0:  points  along  LMX  are  treated 
as  field  points,  with  values  of  the  potential  on  LMX  + 1 set 
equal  to  their  values  on  L = 2.  Thus,  Equation  (30)  is  written  as 


W *!.«>«  N*l. 

<Pk 


n:  /.hx-i 


LHX  u ± 

<t>t 


comparison  with  Eq  (31)  then  leads  to  the  explicit  expression: 


= SYUMX) 


-l-HX  u,3  tux 

Du  - R^  <p^  - C„  FFClmk-}) 


fl*"'”'  i- 


ct""*  EE  iLMx-0 


(b)  Periodicity  condition,  for  E>  here  the  same  formulas 

are  used,  except  that  ig  replaced  by  , 

and  there  is  a contribution  to  from  the  radial  derivative 
of  the  circulation: 
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4> 


LMy( 


= - ^ 


A(t>CN+  1)  - ZA4>(N)  + A(i>CN-i) 

(Ap)^ 


a4>(N+i)  - A(l>(N-i) 

^p^p 


J 


- a4>(n) 


(43) 


It  should  be  noted  that  the  evaluation  of  the  mixed  derivative  at 
LMX,  and  at  the  stations  immediately  off  the  blades  (K  = KLEP 
and  KTEO)  uses  differences  across  the  blade  surface. 

(c)  Blade-surface  condition  for  IBC  = 1:  the  derivative  is 

approximated  by 


^Cr,) 


-i 

A)^;  [ 


<Pk  ^ g <Pk  - 7 <i>K 
iz  Az; 


* i i.l 
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A? 


~ (P^  -h  3 (p^  - 7 <P^ 

IZAK 


+ E. 


^ dsEn 
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(44) 


Using  this  expression,  the  potential  equation  at  LMX  is  written 


as: 


- 2 
T ^ 


+ K-,  Mt"-”  ■[  f,.,  (z  "0 “« - ^ J (4...  _ ^1 
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4 At; 
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This  is  written  in  the  form 


fl£ + c£  = 
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and  is  solved,  along  with  Equation  (31),  to  give 


= svaAix)  = 


CE  -h  BE  ■ FF(lmx-i) 


RE  + BE  . EE(LHX-I) 

(d)  Blade-surface  condition  for  IBC  = 2:  here  the  derivative  cf> 


(47) 


is  approximated  by  using 


P"  AC, 


A ) = A ) ^ - pi  CS) 

tMX  'l=i  ^+P  ^ ' 


Ft; 


(48) 


This  equation,  substituted  into  the  first  form  of  Equation  (44), 
leads  to  the  following  form  of  the  potential  equation 
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The  loading  and  thickness  values  are  calculated  in  BVAL,  and 
are  stored  in: 


0N05  ( K-  KLEP,  N^z) 


-BL 

1 +p^ 


fit  (S') 


(50) 


This  line  relaxation  procedure  is  carried  out  up  to  KMX  - 1.  Then 
values  of  the  potential  at  KMX  are  set  equal  to  the  value  required  by  mass 
conservation: 


C-DRF  ■ - OSf  • V, 


' KHX - 1 


KMX 


DCF 


(51) 


where  the  following  coefficients  are  calculated  in  GRID 


OflF 

DBF 

DCF 


^ ^ KMX  ~ ^ KMX  -l')  / ^ ^ KMX-1  ~ ^ KMX  -z")  KMX  ~ -z) 

C^KMX.  ~ ^ KMX-Z^  ^ ^ Kmx-1  ~ ^ KMX-z)  KMX -1  ~ ^ KMx) 

C^Z  kMX  ~ ^ KMX  - 1 ~ 2 KMX-X  ) / C^KMX  ~ ZkmX  -z)  KMX  ~ ^KMX-1  ) 

(52) 


Next,  the  circulation  is  updated,  using 

A%(^)  = OBV  A<i>*  -h  (l -OBV)  A<I>(N) 

where 

A<f>*  * “irC^rE  ~ 

(Z^-Z^)(ZZ,^-Z^-Z,) 

Zre  = 

^4>a..b  - P yO)  - i(^o..b  > Py^) 


(53) 


(54) 


This  completes  a sweep  in  the  Z-  direction.  The  iteration  counter  ITK  is  then 
incremented,  and  compared  with  ITKMX,  the  maximum  number  of  Z - sweeps  for 
each  radial  one.  After  ITKMX  sweeps  are  done,  the  solution  moves  to  the  next 
radial  station,  where  the  process  is  repeated  for  N = 2,  3,  ...  NMX  - 1. 

The  iteration  counter  ITR  is  used  to  number  the  radial  sweeps;  a total  of 
ITRMX  is  allowed. 


After  each  radial  sweep,  values  of  ^ and  A(p  are  updated  at  N = 1 
and  N = NMX  according  to 


4 

3 


1_ 

3 


fJMA,  L 

Cp 

K 


4 A/MX  - 1 L 

1 


3 


A/MX -a 


€ 


(5S) 


A<p  (i)  = ZA<I>CZ)  - A^(S) 

/:\4j(a/MX)  = ZA^iNMK-t)  - L(^CNMX-Z)  (56) 

Also  updated  at  the  end  of  each  radial  sweep  is  the  quantity  C,  used  [see 
Equation  (10)]  in  enforcing  mass  conservation  at  the  downstream  edge  of  the 
grid.  This  quantity,  called  CONST,  is  evaluated  by  the  trapezoidal  rule. 

The  OUTPUT  subroutine  can  be  called  at  intervals  of  JPRT  in  the  ITK 
iterations  and  at  intervals  of  NPRT  in  the  ITR  iterations.  The  former  choice 
is  intended  for  diagnostic  purposes.  It  would  normally  be  used  only  for 
checkout  purposes,  and  even  then  for  very  few  iterations  and  very  few  radial 
stations,  since  it  generates  many  lines  of  output.  To  avoid  this  diagnostic 
output,  set  JPRT  > ITKMX. 
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Several  options  are  available  in  the  output;  these  are  controlled  by 
the  indicator  lOP,  and  are  described  in  Appendix  A.  The  usual  sequence  is  to 
run  the  solution  for  a certain  number  of  iterations,  write  the  solution  on 
tape  (ISAVE  = 1),  examine  the  output,  and  then  restart  the  solution  (ISTART  = 1) 
with  adjusted  values  of  the  relaxation  factors,  iteration-count  parameters,  etc. 
Some  of  the  output  options  allow  a minimum  number  of  lines  to  be  printed 
during  these  intermediate  stages. 

One  special  option  is  lOP  = 4,  which  calculates  Mach  number  contours. 
This  option  can  only  be  used  with  a solution  found  on  a previous  run:  no 
iterations  are  made,  and  the  Mach  number  calculation  is  done  on  the  values 
read  in  from  the  tape. 

The  quantities  calculated  and  displayed  by  subroutine  OUTPUT  include 
the  perturbation  velocities  and  the  local  Mach 

number.  These  are  defined  as  follows; 


^ - K.j 

^ ^ Cl  + 2Ar 


V 


K 


A/-f 


0 


u 

K 


2Ap 


(57) 


(58) 


j_  _ p dd 

P 8 1;  1+p^  di 


A/  ,£.+  ( A/  . t.  A/,t  NjL-t 


(59) 
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The  first  term  in  this  expression  is  equal  to  d^/d^  , with  a truncation 
error  OCAl^)^,  This  particular  form  is  shown  so  that; 


u- 


r> 


Vi  U-1 

<p 


K + 1 


r; 


2 'Vl+p^' Ai 


(60) 


when 


pA^  _ 7+/0^ 

AZ  ~ p 


(61) 


An  alternate  expression  for  d4>/3i^  . having  the  same  truncation  error,  was 
used  in  References  1-4,  namely; 


-L  It  . 

P dK  zpAn; 


(62) 


This  formula  tends  to  give  erratic  results  when  shock  waves  are  present,  since 
it  maximizes  the  amount  of  differencing  done  across  the  shock. 


The  local  Mach  number  is  defined  by  equating  the  coefficient  of  <p. 


ss 


in  Equation  (4)  to  1-M  : 

i 


, - AJ  = ) - (i  +p  ) - + i) 


or 

Thus 

M VTHi^+O  a,/w;  * ^ ( 

The  square-root  formula  was  used  in  References  1-4;  the  present  report  uses 
the  linearized  version. 
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At  the  very  end  of  every  calculation,  the  output  subroutine  is  called 
with  lOP  = 5.  This  displays  the  following  performance  data  at  each  radius: 

the  values  of  N,  pfN),  A<f)(N) 


vJ/a  ) 

aO  / g 


-5 

2np 


A4> 


(64) 


(65) 


the  turning  angle,  defined  as  ta^ 


-1 


- P — 1 


•oJZ 


1 +/ 


(66) 


the  total  pressure  ratio  'Hz 

'Pot 


= / + 


B 


+ ^ (67) 


2. 


A <f> 


The  convergence  of  the  solution  can  be  monitored  by  calculating  the 
residuals.  This  is  done,  at  ITR  intervals  of  IRXP,  in  subroutine  RESID. 

This  subroutine  evaluates  all  five  terms  in  the  potential  equation  at  the 
interior  points  of  the  grid,  i.e.,  for  N = 2,  3,  ....  MMX  - 1,  K = 2,  3, 

KMX  - 1,  and  L = 2,  3,  LMX  - 1.  The  five  terms  (cal led @ through are 

evaluated  as  follows  [compare  with  Equation  (4)]; 
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(69) 


(At) 


© 


-20  +p'‘)CAr) 


X X 

Mr/  4, 

P f 


tfzr 


J_ 

dp 


(P4>p) 


The  residual  at  each  point  is  defined  as 


L 


@ + (Dr©<-@+® 


(70) 


(71) 


(72) 


(73) 


At  each  value  of  N,  the  value  of  the  maximum  residual  is  printed,  along  with 
the  K and  L values  of  the  point  where  it  occurs.  Also  shown  is  the  average 

V L 

residual,  i.e.,  the  mean  of  over  the  range  N = 2,  NMX  - 1,  K = 2, 

KMX  - 1;  L = 2,  LMX  - 1.  The  magnitudes  of  the  residuals  themselves  are  only 
meaningful  in  relation  to  the  sizes  of  the  individual  terms  whose  sum  they 
represent.  Thus  the  mean  value  of  <p  and  of  the  sum  of  the  absolute  values  of 
the  five  terms  are  also  shown,  i.e.: 


PHI  = ^ z:  I V 

/V  =-2  , - f 

K =2,  KMX.-1 
L =2  , X - / 


(74) 


20 


SUM  . 1 2.^  {|®|*|®I*I©I*I®M®||  (7s, 

K = a,  KMx-1 

X.  = a,  LMx  -I 

where  P = (NMX  - 2)  (KMX  - 2)  (LMX  - 2)  (76) 

A second  way  to  monitor  the  convergence  of  the  solution  is  to  observe 
the  perturbation  velocities  at  a number  of  selected  points.  Our  experience 
has  been  that  these  will  usually  converge  to  three  decimal  places  when  AVG  RES/ 
AVG  SUM  C 3 X 10‘^. 


If  IRXP  is  prefixed  by  a minus  sign,  subroutine  RESID  will  display 
the  values  of  all  five  terms  and  the  residual,  at  each  of  the  interior  grid 
points.  This  option,  which  produces  many  lines  of  output,  can  be  used  for 
diagnostic  purposes. 


Section  3 
SAMPLE  CASE 


To  illustrate  the  operation  of  the  program,  some  details  are  presented 
of  the  two-dimensional  case  shown  in  Figure  7 of  Reference  4.  This  cascade  had 
a stagger  angle  of  45°,  an  inlet  relative  Mach  number  of  0.9,  a solidity  ^a./^r 
of  0.411,  six  percent  thickness-to-chord  ratio,  and  no  camber  or  thickness.  The 
input  data  were  as  follows: 

B = to  , A ^ , C./Lj.  = 0.<lo73  , = o.  , M = O-U’fa’li 

i = 1.4  , FXB  = -Z  , FXI  = +z  , PXE  = 7 , RXH  = 0.9 

OBV  = 0.05 

KMX  ^ UO  , LMX  - 30  , NMX  = 3 , IBC  =■  l , IDM  = 2 , IT  KMX  = 3L,o  , 
JPRl  = 400  . ISTRRT  - O . ISftVE  - O , lOP  = / , IT R P\X  = i , 

NPRT  ^ 0 , I TPR  = / , IRXP  ^ 1 

5/(7)  = 0.089S5Z3  , 5/(2;  throu<^h  BVOo)  = <0- 


Figures  3a  - 3c  arc  the  first  three  pages  of  the  output,  showing  the 
run  conditions,  coordinate  system,  and  blade  geometry.  Figures  4a  - 4d  are 
portions  of  the  pages  on  which  the  potential,  streamwise,  normal  and  radial 
velocity  components  appear.  The  last  of  these  is  zero  for  the  case  shown;  it 
would  be  non-zero  if  NMX  were  greater  than  3,  i.e.,  if  two-dimensional  strip- 
theory  calculations  were  being  done  at  a series  of  radial  stations.  Figure  5 is 
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the  last  page  of  the  output,  showing  the  most  recent  value  of  the  circulation, 
the  residuals,  and  the  performance  data.  Note  that  the  average  residual  is 
already  down  to  less  than  1/100  of  the  average  sum  of  terms,  even  at  this  early 
stage  of  the  iterations. 

This  case  used  3.5  minutes  of  CPU  time  on  an  IBM  360/65. 
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Figure  1 BLADE-FIXED  COORDINATES 

ROTOR  IS  STATIONARY  IN  A 
HELICAL  APPROACH  FLOW 


Figure  2 DEFINITIONS  OF  BLADE-SURFACE  GEOMETRY 
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SAMPLE  CASE  FOR  USERS  GUIDE  REPORT 

CASE  IS  THAT  OF  FIGURE  7.  JOUNAL  OF  ENERGY,  P.  289,  SEPT/OCT  1977 
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FIGURE  3a 


SAMPLE  CASE  FOP  USERS  GUIDE  REPORT 
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FIGURE  3a  (Cont. 
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FIGURE  3b 
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VALUES  OF  THE  POTENTIAL  AFTER  ITR  - 1,  ITK  • 360  AT  RHO(  2)  - I. 0000 
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FIGURE  4a 


VALUES  OF  THE  POTENTIAL  AFTER  ITR  « 1,  ITK  • 360  AT  RHO(  2)  - 1.0000  (Cont . 
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FIGURE  4a  (Cont. 
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FIGURE  4c  (Cont. 


NOTE:  THE  FOLLOWING  VALUES  OF  V RADIAL  ARE  BASED  ON  THE  2D  STRIP-THEORY  APPROXIMATION 
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FIGURE  4d 


MOTE:  THE  FOLLOWING  VALUES  OF  V RADIAL  APE  BASED  ON  THE  20  STRIP-THEORY  APPROXIMATION  (Cont. 
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FIGURE  4d  (cont. 
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APPENDIX  A 

GUIDE  TO  PREPARING  THE  INPUT 


Cards  I and  2:  FORMAT  I8A4:  These  cards  are  for  run  identification;  any 

alphameric  characters  may  be  used.  Neither  card  can  be  omitted. 

Card  FORMAT  5E15: 


Columns  1-15;  B,  the  number  of  blades.  For  a two-dimensional  cascade 
solution,  use  B=10  (and  see  comments  below  regarding  1UM=2) . 

Columns  16-30:  h,  the  hub-to-tip  radius  ratio.  For  a two-dimensional 
cascade  solution,  use  h = 0.9. 

Columns  31-35:  Solidity,  C^/ i-j  = / — g — ' 

Columns  4fa-60:  Axial  Mach  number  at  the  inlet,  or*  / O-no 


Columns  61-75:  Tangential  Mach  number  at  the  tip,  OOC^/o.^.  The  present 
version  of  the  program  was  developed  for  inlet  relative  Mach  numbers 
/ o-go  less  than  one.  The  program  cannot  handle  values 


greater  than  one,  because  the  boundary  conditions  imposed  at  the  upstream 
and  downstream  edges  of  the  grid  are  not  radiation  conditions. 


Card  4:  FORMAT  5E15 


Columns  1-15:  ,the  specific-heat  ratio. 

Columns  10-30:  FXB:  The  upstream  edge  of  the  grid  will  be  located  at  FXB 
times  the  axial  projection  of  the  chord  . Values  of  FXB  = -2  or  -3 

have  been  used.  FXB  must  be  negative. 

Columns  31-45:  FXl : The  downstream  edge  of  the  grid  will  be  located  at 
FXl  • downstream  of  the  trailing  edge.  FXI  = -FXB  will  place  the 

greatest  concentration  of  grid  points  near  the  mid-chord  of  the  blades. 

Columns  46-60:  RXE  is  the  relaxation  factor  used  at  subsonic  points. 

Values  up  to  1.2  (and  "tapered",  i.e. , ITPR=1  - see  below)  have  been 
used  for  the  off-design  problem  (IBC=1),  but  very  little  is  known  about 
the  limitations  on  this  parameter.  For  the  design  problem  (1BC=2) , it 
may  be  necessary  to  use  values  as  low  as  0.1,  with  1TPR=0,  early  in  the 


A-1 


iterations,  increasing  to  0.8  (and  rn’R=0)  at  the  later  stages.  Again, 
very  little  is  known  about  these  limitations. 

Columns  61-75:  RXH  is  the  relaxation  factor  used  at  supersonic  points. 
RXH=0.9  has  been  used  successfully.  The  ITPR  parameter  does  not  affect 
RXH  - a constant  value  is  used,  at  all  supersonic  points  in  the  grid. 

Card  5:  FORMAT  5E15 

Columns  1-15:  OBV,  the  relaxation  factor  used  in  updating  the  circulation. 
Values  of  0.05  and  0.1  have  been  used,  but  the  limits  on  this  parameter 
and  their  coupling  to  the  other  relaxation  factors  are  unknown. 

Card  6:  FORMAT  1615 

Columns  1-5/6-10/11-15:  Grid  size  parameters  KMX/LMX/NMX.  These  are 

limited  to  60/30/10,  by  a COMMON  statement  in  the  version  listed  below. 
The  number  of  calculations  done  in  a given  S--sweep  is  proportional  to 
KMX,  and  the  number  of  Z-sweeps  required  for  convergence  is  also 
proportional  to  KMX  (since  it  takes  KMX  iterations  before  information 
at  the  downstream  edge  of  the  grid  affects  the  solution  at  the  upstream 
edge).  Thus,  the  total  time  required  for  a given  calculation  is  likely 
to  vary  as  the  square  of  KMX,  and  values  as  low  as  possible  should  be 
used . 

If  shock  capturing  is  important,  LMX  should  then  be  chosen  so  as  to 
make  approximately  equal  to  p at  a radius  in  the 

center  of  the  region  where  shock  definition  is  desired.  For  the  parti- 
cular stretching  of  the  ^-coordinates  used  here,  and  for  the  special 
case  where  FXI  = -FXB,  the  value  of  p^t^  /AE  at  midchord  is 

±T_  J_  (KMX  -I-  1)  JU  KMX 

Ft  ' ' 2.(1  *Z  FXI)  ' LMX  - ; 

This  equation  can  be  solved  for  LMX,  once  the  other  parameters  are 
chosen. 

The  quantity  NMX  must  be  three  or  greater;  for  a two-dimensional 
cascade  solution,  set  it  equal  to  3. 

A- 2 


tJ 


Column  20:  IBC=1  for  the  off-design  case,  IBC=2  for  the  design  case. 
Column  25:  1DM=3  for  a three-dimensional  solution.  For  a two-dimensional 


cascade  solution,  set  IDM=2  and  NMX=3.  For  this  case,  the 
hub/tip  ratio  -k  and  the  number  of  blades  B are  meaningless,  and  can 
be  read  in  as  any  arbitrary  numbers;  the  values  h = 0.9  and  B = 10  arc 
recommended.  The  desired  cascade  parameters  are  actually  the  inlet 
relative  Mach  number  and  the  stagger  angle.  These  determine  p , 
which  is  the  ratio  of  the  pitchwise  to  axial  velocity  components  upstream 
of  the  cascade: 


p = U>r  / a. 

• A 


The  hub  and  tip  radii  are  then  determined  by  the  fact  that,  for  NMX=3, 
the  value  of  p lies  half-way  between  /’h  and  : 


P = 


- £ilL£i  - 


Pr  i 


Thus, 


Z 2-k. 

^ / + A P ^ ~ 77X  P 


{k-l^ 


The  desired  value  of  the  pitchwise  component  of  the  inlet  Mach  number 
at  station  p is 


= . P M ^ 

Pr 


Pr 


M = 


P 


(A- 3) 


+ ^ " 1+A  V/  +p 

For  a two-dimensional  strip-theory  calculation  of  a three-dimensional 
case,  set  II)M=2  and  3 < MMX  ^ 10.  In  the  latter  case,  all  radial 
derivatives  are  set  equal  to  zero,  and  each  radial  station  is  treated 
as  though  it  were  a two-dimensional  section. 


Columns  26-30:  ITKMX  is  the  maximum  number  of  iterations  that  will  be 
made  in  the  S-direction,  on  each  radial  sweep. 

Columns  31-35:  .JRPT : Output  is  printed  at  intervals  of  JRPT  in  the  •2'- 
direction  iterations  on  every  radial  sweep.  It  is  used  only  for  looking 


A- 3 


at  intermediate  results,  and  gives  a large  amount  of  output.  Normally, 
results  are  needed  only  after  a given  number  of  radial  iterations  (see 
NRPT  below).  To  suppress  these  intermediate  results,  set  .JRPT  ? ITKMX. 

Columns  36-40:  ISTART  = 0 if  you  are  starting  from  scratch;  ISTART  = 1 if 
you  are  reading  in  values  of  <p  (on  tape)  from  a previous  run. 

Columns  41-45:  1SAVE=1  will  write  all  values  of  A and  cf>(L,K^N)  on 

tape  at  the  end  of  the  calculation.  1SAVE=0  writes  no  tape. 

Columns  46-50:  lOP  controls  what  is  printed  on  output.  See  the  comment 
cards  in  subroutine  OUTPUT  for  details. 

Columns  51-55;  ITRMX  gives  the  maximum  number  of  radial  iterations.  On 
each  radial  iteration,  ITKMX  axial  iterations  are  done,  at  radial 
stations  N=2,  NMX-1. 

Columns  56-60:  NRPT:  Output  will  occur  at  intervals  of  NRPT  in  the  radial 
iterat ions . 


Columns  61-65;  1TPR=1  will  cause  the  relaxation  factor  used  at  elliptic 
points  to  be  "tapered",  according  to 


6<J 


1 + (RXB  - 1) 


a 


cue 


(A-4) 


1TPR=0  usesto=RXE  at  all  elliptic  points. 

Columns  66-70:  Residuals  are  calculated  and  printed  at  intervals  of  IRXP 
in  the  radial  iterations.  If  IRXP  is  preceded  by  a negative  sign,  the 
parameter  ISHO  is  internally  set  equal  to  1,  and  absolute  values  of  all 

terms  in  the  potential  equation  are  displayed  as  well. 

Cards  7 and  8:  These  contain  the  input  data  for  blade  shape  and  loading  - see 

comments  in  subroutine  BVAL  for  details.  Both  cards  must  be  used, 
with  dummy  zeros,  if  necessary. 

Card  8:  If  I0P=4,  this  card  specifies  the  locations  and  Mach-number  values  at 


which  Mach-number  contours  are  calculated. 


APPENDIX  B 
PROGRAM  LISTING 


LEVEL  21.7  < DEC  72  ) 


OS/360  FORTRAN  H 


ISN 


ISN 

ISN 

ISN 

ISN 

ISN 


ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 


COMPILER  OPTIONS  - NAME=  MAI N . OPT=02 , L 1 NECNT=60 , S IZE=0000K , 

SOURCE , EBCDIC, NOLI  ST, NOOECK, LOAD , MAP ,NOEO IT, ID, XREF 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM  3DNSDA  - THREE-DIMENSIONAL  NONLINEAR  SMALL-D I STURBANCE 
(VERSION  A)  PROGRAM  FOR  CALCULATING  THE  FLOW  THROUGH  AN 
ISOLATED  COMPRESSOR  BLADE  ROW.  DOCUMENTATION  IS  GIVEN  IN: 

W.  J.  RAE,  NONLINEAR  SMALL -D I STURBANCE  EQUATIONS  FOR 
THREE-DIMENSIONAL  TRANSONIC  FLOW  THROUGH  A 
COMPRESSOR  BLADE  ROW,  CALSPAN  CORPORATION 
REPORT  AB-5487-A-1,  AFOSR  TR-76-1082, 

AD  A031234  (AUG  1976) 

- - - RELAXATION  SOLUTIONS  FOR  THREE-DIMENSIONAL 

TRANSONIC  FLOW  THROUGH  A COMPRESSOR  BLADE  ROW 
IN  THE  NONLINEAR  SMAL L -D I STURBANCE  APPROXIMA- 
TION, CALSPAN  CORPORATION  REPORT  AB-5487-A-2, 
AFOSR  TR-76-1081,  AD  A032B53  (AUG  1976) 

- - - CALCULATIONS  OF  THREE-DIMENSIONAL  TRANSONIC 

COMPRESSOR  FLOWFIELDS  BV  A RELAXATION  METHOD, 
AIAA  JOURNAL  OF  ENERGY,  VOL.  I,  NO.  5. 

SEPT  - OCT  1977,  PP  284  - 296 

USERS  ARE  URGED  TO  COMMUNICATE  WITH  THE  AUTHOR,  CONCERNING 

PROBLEMS,  RESULTS,  AND  SUGGESTED  CHANGES  - 

AERODYNAMIC  RESEARCH  DEPARTMENT,  CALSPAN  CORPORATION 

PO  BOX  235 

BUFFALO,  NY  14221 

TELEPHONE  716/632-7500 


0002 


0003 

0004 

0005 

0006 
0007 


0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 
0019 


COMMON  ZT( 30 ),XI(60 ),RO{ 10  ),F(  30,60, 10  ),A( 10  ) ,B( 10)  ,D( 10),E( 10  ), 

• EE(30),FF(30),SV(30),DPHI(10>,BV(10),US(30,60),UN(30,&0), 

• DNDS( 30, 10,2  ) ,FAK( 60  ) ,FHW( oO  ) ,UTA( 30  ) ,UTB( 30  ) ,UTC( 30  ) ,BBV( 30  ) 

• BSA{ 30 ) ,BSB( 30 ) ,BSC( 30  ) ,VSA( 30  ) , VSB( 30  ) ,BAV( 30  ) ,LOP( 30  ) 

• ,SONIC(10) 

COMMON  DAF ,DBF ,DCF , AKLKE , RX , OX , CONST 

COMMON  H , DRO , DZT , TDZ ,TPP,XIC,XIB,XID , TDR , AKLK , DR02 , 0ZT2 
COMMON  OMCU,RHUB,RTIP,R2,RXE ,RXH,UA,UB,TRDZ,ERA,ERB,CKLKE,SB, 

• BTR,FXB,FXI , EM2 , AMA , AMB , DM , DKA , DKR , DO , DOS , TZ , TTZ , TZS , TRS , TDQ 
COMMON  K,N, IDM,  KMX , LMX , NMX , KMXM 1 ,KMXM2,LMXM1 ,LMXM2,NMXM1 . 

• NMXM2,KTE0, IRX,  I TK , ITR , I TRMX , KLE P , I TPR , KUP , KDN 
COMMON  ID(36) 

C 

C THE  FOLLOWING  CALSPAN  ROUTINE  PLACES  A ZERO  IN  THE  LOCATIONS 

C FROM  ZT( 1 ) THROUGH  ID(36) 

C 

CALL  CLEAR(ZT( 1 ) , ID( 36  ) ) 

C 

READ! 5, 102  ) ( I D( I ) , I = 1 , 36  ) 

102  FORMAT! 18A4  ) 

WRITE!  6,200  ) 

200  FORMAT!  IHl  ) 

WRITE! 6,201  ) ! I D ! I ) , ! = 1 , 36  ) 

201  FORMAT! 30X, 1PA4  ) 

READ! 5, 100  ) BN , H , CALT , EMX , EMTG , GMA , FXB , FX I , RXE , RXH ,OBV 

100  FORMAT! 5E15. 5 ) 

READ! 5, 101  ) KMX, LMX, NMX, IBC, IDM, ITKMX, JPRT, I START, I SAVE , I OP 

• , ITKMX, NPRT, ITPR, IRXP 

101  FORMAT! 1615) 

READ! 5, 100  ) ! BV! I ) . I = 1 , 1 0 ) 


B-1 


ISN 

0020 

IBN  = INK  BN  + 0.01  ) 

ISN 

0021 

WRITE(6.202>  IBN.H.CALT 

ISN 

0022 

202  FORMAT! //5X. ‘THIS  BLADE  ROW  HAS', 14,’  BLADES,  WITH  A HUB-TO-TIP' 
* ‘ RATIO  OF  ',F5.3.'  AND  SOLIDITY  CA/LT  = ',F5.3) 

ISN 

0023 

RTIP  = EMTG/EMX 

ISN 

0024 

RHUB  = H*RTIP 

ISN 

0025 

EM2  = EMX'EMX 

ISN 

0026 

REi.  = SORT!  EM2  + EMTG»EMTG  ) 

ISN 

0027 

WRITE(6,203)  EMX , EMTG , REL ,GMA 

ISN 

0028 

203  FORMAT!  /5X. 'AXIAL  MACH  NO.  » ',F5.3,',  TANGENTIAL  MACH  NO  AT' 

• ' THE  TIP  = ',F5.3,',  TOTAL  MACH  NO.  AT  THE  TIP  = ',F5.3, 

• //5X, 'SPECIFIC  HEAT  RATIO  =‘.F6.3) 

ISN 

0029 

IF! REL  .GT. 1 . > WRITE!  6 ,207  ) 

ISN 

0031 

207  FORMAT! //15X. 'WARNING  - INLET  RELATIVE  MACH  NUMBER  AT  THE  TIP  ', 

* 'IS  SUPERSONIC . 

* /15X. ■ BOUNDARY  COMDITIONS  UPSTREAM  AND  DOWNSTREAM  DO  NOT  ', 

“ 'CONTAIN  THE  REQUIRED  RADIATION  CONDITIONS', 

* /15X, 'SOLUTION  UPSTRLAM  MAY  BE  INVALID', 

* / 15X ,' SOLUTION  DOWNSTREAM  MAY  NOT  CONVERGE,  ESPECIALLY  IF', 

* ' THE  OUTLET  RELATIVE  MACH  NUMBER  IS  SUPERSONIC',/) 

ISN 

0032 

C 

C 

IF!  ISTART.NE. 1 ) GO  TO  98 

READ  STARTING  VALUES  FROM  TAPE 

ISN 

0034 

KS  = KMX 

ISN 

0035 

LS  = LMX 

ISN 

0036 

NS  = NMX 

ISM 

0037 

READ! 1 ) LMX. KMX, NMX, DPMI 

ISN 

0038 

DO  97  N=1,NMX 

ISN 

0039 

READ! 1)  !!F!L,K,N),L  = 1,LMX),K=1 ,KMX  ) 

ISN 

0040 

97  CONTINUE 

ISM 

0041 

IF!KS.NE.KMX>  GO  TO  33 

ISN 

0043 

IF!LS.NE.LhX)  GO  TO  33 

ISN 

0045 

IF! NS.NE . NMX  ) GO  TO  33 

ISN 

0047 

GO  TO  93 

ISN 

0048 

33  WRITE!6,210)  KS , LS , NS , KMX , LMX , NMX 

ISN 

0049 

C 

c 

210  FORMAT! //5X, 'WARMING  - GRID  SIZE  ON  CARD  INPUT  DIFFERS  FROM', 

* ' THAT  ON  TAPE ' , 

* //10X,'CARD  KMX  = ',I4,‘  LMX  = ',14,'  NMX  = ',14, 

* //I  OX,  'TAPE  KMX  ^ ',14,'  LMX  = ',14,'  NMX  = ',14,) 

END  OF  TAPE  INPUT 

ISN 

0050 

L 

98  WRITE!  6,204  ) I DM , KMX , LMX , NMX 

ISN 

0051 

204  FORMAT!  /5X,'THIS  IS  A ' , 1 2 , ' -0 1 ME  NS lONAL  CALCULATION,  WITH  GRID' 
• ' SIZE  KMX/LMX/NMX  = ’ , I 3 , ' / ' , I 3 , ' / ' , I 3 ) 

ISN 

0052 

WRITE! 6 ,206 ) 

ISN 

0053 

206  FORMAT!/5X, 'RELAXATION  FACTORS  FOR  ELLIPTIC  AND  HYPERBOLIC, 

* ' POINTS  ARE  LISTED  BELOW  AS  RXE  AND  RXH’> 

ISN 

0054 

SB  = 1 GMA  + 1 . >«EM2 

ISN 

0055 

TPI  = 2.0*3.14159265 

ISN 

0056 

IPB  = TPI/BN 

ISN 

0057 

OMCU  = T!’a*RTIP*CALT 

ISN 

0058 

LMXMl  = LMX  -1 

ISN 

0059 

LMXM2  = 1 MX  - 2 

ISN 

0060 

KMXMl  = KMX  - 1 

ISN 

0061 

KMXM2  = KMX  -2 

ISN 

0062 

NMXMl  = NMX  - 1 

ISN 

0063 

NMXM2  = NMX  - 2 

ISM 

0064 

CALL  GRID 

ISN 

0065 

DO  99  I = I, LMX 

ISN 

0066 

99  LOP! I ) = I 

ISN 

0067 

IF!10P.NE.4)  GO  TO  15 

ISN 

0069 

REA0!5,104)  AHA , AMB . DM . KUP , KDN 

ISN 

0070 

104  FORMAT! 3E15. 5, 315) 
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ISN 

0071 

DO  21  N = 1 .NMX 

ISN 

0072 

UA  = RO(N) 

ISN 

0073 

UB  =( 1 . /< 1 .+RO< N >*RO< N ) ) )/TDQ 

ISN 

0074 

TRDZ  = TDZ*RO(N) 

ISN 

0075 

CALL  OUTPUT! lOP . IBC  > 

ISN 

0076 

21 

CONTINUE 

ISN 

0077 

GO  TO  30 

isri 

0078 

L 

15 

CALL  BVAL( IBC ) 

ISN 

0079 

9 

OBM  - 1 . - OBV 

ISN 

0080 

ITR  = 1 

ISM 

0081 

NPT  = 0 

ISN 

0082 

IRXPT  = 0 

ISN 

0083 

ISHO  ^ 0 

ISN 

0084 

IF< IRXP.GE.O)  GO  TO  44 

ISN 

008G 

IRXP  = -IRXP 

ISN 

0067 

ISHO  = 1 

ISN 

0088 

44 

BGL  = RO( 1 )*DPHI(  1 )/2. 

ISN 

0089 

DO  45  N = 2,NMXM1 

ISN 

0090 

45 

BGL  = BGL  + RO(N>*DPHHN) 

ISN 

0091 

BGL  = (BGL  + RO( NMX >»DPHI( NMX >/2 

ISN 

0092 

CONST  =-2.*BGL/(TPB*(  1 .-EM2  )»<RT 

c 

THE  FOLLOWING  STATFMENTS  ARE  USED 

c 

c 

WHEN  THE  SOLUTION  IS  BEING  FOUND  A 

ISN 

0093 

NA  = 1 

ISN 

0094 

NB  = NMX 

ISN 

0095 

IF(NMX.NE.3>  GO  TO  304 

ISN 

0097 

NA  = 2 

ISN 

0093 

NB  = 2 

c 

c 

BEGINNING  OF  RHC  - SWEEP 

ISN 

0099 

304 

DO  300  N = 2,NMXM1 

ISN 

0100 

R2  = RO(N)»RC(N) 

ISN 

0101 

BTR  =1.0-  EM2*( 1 .0+R2 ) 

ISN 

0102 

AKLK  = D(N)*TZS 

ISM 

0103 

UA  = RO(N) 

ISN 

0104 

UB  =(  1 ./(  1 . + RO(N>*RO(N>n/TDQ 

ISN 

0105 

TRDZ  = TDZ*RO(N) 

ISN 

0106 

DKA  = -B(N)*TTZ 

ISN 

0107 

ERA  = 1 . + 0 . 5*DRO/RO( N ) 

ISN 

0108 

ERB  = 1.  - 0 . 5*DRO/RO( N ) 

ISN 

0109 

ITK  = 1 

ISN 

0110 

IPRT  = 0 

ISN 

0111 

IF( IDM.NE .2  ) GO  TO  1 

ISN 

0113 

CONST  =-DPHl( N )/(TPB»< 1 .-EM2 ) ) 

L 

C 

Q 

BEGINNING  OF  XI  - SWEEP 

ISN 

0114 

I 

DO  2 L » 1 ,LMX 

ISN 

0115 

UTA(t.  > = 0.0 

ISN 

0116 

UTB<  L ) = 0.0 

ISN 

0117 

VSA(L)  = 0. 

ISN 

0118 

BAV( L ) = BTR 

ISN 

0119 

2 

CONTINUE 

ISN 

0120 

38 

DO  5 K = 2,KMXM1 

ISN 

0121 

RX  = I.  + (RXE-1 . )»EXP(-{<XI(K)-( 

ISN 

0122 

IF( ITPR.EQ.O  ) RX  = RXE 

ISN 

0124 

OX  = I.  - l./RX 

ISN 

0125 

RXM  = 1.0  - RX 

ISN 

0126 

IRX  • 1 
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ISN 

0127 

C 

C 

C 

BEGINNING  OF  ZETA  - SWEEP 

4 AKLKE  = AKLK/FAK<K> 

ISN 

0128 

DKR  = TRS*E( N )/FAK( K ) 

ISN 

0129 

CKLKE  = AKLKE 

ISN 

0130 

29  DO  3 L = 1 ,LMX 

ISN 

0131 

UTC( L > = UTB<  L ) 

ISN 

0132 

UTB( L ) = UTA<  L > 

ISN 

0133 

UTA<  L > = F( L ,K.N  ) 

ISN 

0134 

BBV( L ) = 3AV( L > 

ISN 

0135 

BAV(L)  = BTR  - SB*FAK( K )*{F{ L ,K+1 .N  }-F< L,K-1 ,N  ) )/TDQ 

ISN 

0136 

VSB(  L > = VSA<  1.  ) 

ISN 

0137 

VSA(L>  = 0. 

ISN 

0138 

IF( BAV( L > .GT.O . > GO  TO  8 

ISN 

0140 

VSA(L)  = 1. 

ISN 

0141 

IRX  = 2 

ISN 

0142 

8 8SA(L>  = ( 1 .-VSA( L ) )*{ BAV{ L >+A< N ) > 

ISN 

0143 

BSB( L > = VSB{ L )*BBV( L > 

ISN 

0144 

BSC( L > = VSA( L )«A( N ) 

ISN 

0145 

r- 

3 CONTINUE 

ISN 

0146 

L 

c 

c 

c 

UTA.UTB,  AND  UTC  NOW  CONTAIN  DATA,  FROM  THE  PREVIOUS  ITERATION, 
FOR  F( L .K.N > ,F( L ,K-1 .N  ) , AND  F(L,K-2,N>.  RESPECTIVELY 

14  CALL  BVAL0<  IBC  > 

ISN 

0147 

DO  7 L = 2,LMXM1 

ISN 

0148 

BKL  = -2. "AKLKE  - 2 . "FAKI K ) "BSAI L ) /RX  + 2 . "F HW< K- 1 >*BSB( L > 

ISN 

0149 

* - BSC( L )*FHW( K-1  ) 

DKL  = -BSA( L )*< FHW( K )*F( L ,K+1 ,N >-2. »FAK( K )"OX"UTA< L > 

ISN 

0150 

* + FHW(  K-1  )*F(  L .K-1  .N  ) ) i-  BSB<  L >•  ( FHW(  K- !>  * ( UTA(  L ) + F < L . K- 1 

* + FHW{ K-2  )"( F( L ,K-1 .N  )-UTC( L > > ) - BSC  1 L ) * ( FHW( K )•( F ( L , K+ 1 

* -UTa(L>>  + FHW( K-1 >"F( L .K-1 ,N > ) 

* + 2. *DKA*( -F( L + 1 ,K-1 ,N > + F ( L ,K-1 ,N  ) + F( L + 1 ,K,N >-2 .»F< L .K.N ) 

* + F( L-1  .K.N > + F( L ,K+1  ,N  )-F( L-1 ,K+1 ,N  ) ) 

* -DKR"( EKA*F (L,K,N+1)  + ERB*F(L,K,N-1  )-2. "UTAI L ) ) 

DNM  = BKL  + CKLKE*EE( L-1  ) 

ISN 

0151 

EE<L)  = -AKLKE/DNM 

ISN 

0152 

16  FF(L>  = (DKL  - C KL KE *F F ( L - 1 > ) /DNM 

ISN 

0153 

7 CONTINUE 

ISN 

0154 

CALL  BVALK  IBC) 

ISN 

0155 

406  DO  19  J = 1 .LMXMl 

ISN 

0156 

L = LMX-J 

ISN 

0157 

SV( L > = EE( L >*SV( L + 1 ) + FF( L > 

ISN 

0158 

19  CONTINUE 

ISN 

0159 

c 

70  IF( IRX.EQ. 1 > GO  TO  6 

ISN 

0161 

11  IF( IRX.EQ. 2)  RX  = RXH 

ISN 

0163 

p 

RXH  = 1.0  - RX 

ISN 

0164 

6 DO  12  L = l.LMX 

ISN 

0165 

p 

12  F(L,K,N>  = RX*SV(L)  + RXM*F(L.K,N) 

ISN 

0166 

p 

5 CONTINUE 

ISN 

0167 

DO  41  L = l.LMX 

ISN 

0168 

41  F<L,KMX,N)  = {C0NST-DAF*F(L,KMXM2,N>-DBF"F«L,KMXM1,N>>/DCF 

ISN 

0169 

p 

307  IF< IBC.EQ.2>  GO  TO  43 

ISN 

0171 

c 

c 

RESET  DPMI 

47  XI  = XI(KTE0-2) 

ISN 

0172 

X2  = XKKTEO-I  ) 

ISN 

0173 

X3  = XID 

ISN 

0174 

CD  - (X2-X1 >"(2."X3-X1-X2) 

r 
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ISN 

0175 

EBO  = -{ X3-X2 >•( X3-X2 >/CD 

ISN 

0176 

FBD  = ( X3-X1  )•( X3-X1  >/CD 

ISN 

0177 

DFV  = EBD*( F( 1 ,KTE0-2.N  ) - F ( LMX , KTEO-2 . N ) > 

+ FB0»< F< 1 ,KTE0-1  ,N > - F < LMX . KTEO- 1 . N ) > 

ISN 

0178 

DPHKN)  = 06V*DFY  + OBM*DPHI(N) 

ISN 

0179 

L 

43 

IPRT  = IPRT  + 1 

ISN 

0180 

IF( IPRT. LT.OPRT  ) GO  TO  42 

ISN 

0182 

IPRT  = 0 

ISN 

0183 

CALL  OUTPUT! lOP , IBC ) 

ISN 

0184 

42 

ITK  = ITK  + 1 

ISN 

0185 

IF(  ITK.GT. ITKMX > GO  TO  300 

ISN 

0187 

GO  TO  1 

ISN 

0188 

300 

CONTINUE 

ISN 

0189 

u 

306 

IF(NMX.EQ.3)  GO  TO  310 

ISN 

0191 

DO  301  K = l.KMX  i 

ISN 

0192 

DO  302  L = 1 ,LMX 

ISN 

0193 

F<L,K,1)  = (4.*F(<.,K,2)-F<L,K,3>)/3. 

ISN 

0194 

F(L,K,NMX)  = (4.*F(L,K.NMXM1  )-F(L.K,NMXM2>)/3. 

ISN 

0195 

302 

CONTINUE 

ISN 

0196 

301 

CONTINUE 

ISN 

0197 

DPHK 1 ) = 2.*DPHI(2)  - 0PHI(3) 

ISN 

0198 

DPHKNMX)  = 2.*DPHI(NMXM1  ) - DPHI<NMXM2> 

ISN 

0199 

GO  TO  305 

ISN 

0200 

310 

DO  311  K = l.KMX 

ISN 

0201 

DO  312  L = 1 ,LMX 

ISN 

0202 

F(L,K,1)  = F(L,K,2) 

ISN 

0203 

F(L.K.3)  = F(L,K,2> 

ISN 

0204 

312 

CONTINUE 

ISN 

02C5 

311 

CONTINUE 

ISN 

0206 

DPHK  1 > = DPHK  2 > 

ISN 

0207 

D?HI<3>  = DPH1<2) 

ISN 

0208 

305 

NPT  = NPT  + 1 

ISN 

0209 

IF<NPT.LT.HPRT)  GO  TO  62 

ISN 

0211 

NPT  = 0 

ISN 

0212 

ITK  = ITK  - 1 

ISN 

0213 

DO  63  N = NA,NB 

ISN 

0214 

UA  = RO<N) 

ISN 

0215 

UB  =(  1 ./( I .+RO{ N )»RO( N ) ) >/TDQ 

ISN 

0216 

TRDZ  = TDZ*RO(N) 

ISN 

0217 

CALL  OUTPUT! lOP , IBC ) 

ISN 

0218 

63 

CONTINUE 

ISN 

0219 

62 

CONTINUE 

ISN 

0220 

IRXPT  = IRXPT  + 1 

ISN 

0221 

IF! IRXPT.NE . IRXP > GO  TO  32 

ISN 

0223 

IRXPT  = 0 

ISN 

0224 

CALL  RESID! ISHO) 

ISN 

0225 

32 

CONTINUE 

ISN 

0226 

IF!  IBC.EQ.2>  GO  TO  34 

ISN 

0228 

BGL  = RO! 1 )*DPHI! 1 >/2. 

ISN 

0229 

DO  46  N = 2,NMXM1 

ISN 

0230 

46 

BGL  = BGL  + RO! N >*DPHI 1 N ) 

ISN 

0231 

BGl  = !BGL  + RO! NMX  )*DPHII NMX )/2. )*DRO 

ISN 

0232 

CONST  =-2.*BGL/!TPB»! I . -EM2  )•! RT I P ‘RTI P-RHUB*RHUB ) ) 

ISN 

0233 

34 

CONTINUE 

ISN 

0234 

ITR  = ITR  + 1 

ISN 

0235 

C 

IF! ITR.LE.ITRMX)  GO  TO  304 
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ISN  0237 


92  IF(  ISAVE  .NE  . 1 ) GO  TO  30 


C 

C 

r- 

SAVE  VALUES 

ON  TAPE 

ISN 

0239 

WRITE! 2 ) 

LMX, KMX. NMX, DPMI 

ISN 

0240 

DO  95  N= 

1 ,NMX 

ISN 

0241 

WRITE! 2 ) 

! ! F ! L , K , N > . L = 1 , LMX ) . K= 1 , KMX  > 

ISN 

0742 

95  CONTINUE 

c 

c 

END  OF  TAPE 

WRITE 

ISN 

0243 

30  CONTINUE 

ISN 

0244 

CALL  OUTPUT! 5, IBC > 

ISM 

0245 

STOP 

ISN 

0246 

END 
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LEVEL  21.7  ( DEC  72  ) 


OS/360  FORTRAN  H 


ISN 

ISN 


ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 


COMPILER  OPTIONS 
C 


NAME=  MAIN,OPT=02.LINECNT=60,SIZE*0000K. 
SOURCE,EBCDIC,NOLIST,NOOECK,LOAO,MAP,NOEDIT,ID,XREF 


0002 

0003 


0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0044 

0045 

0046 

0047 

0048 


SUBROUTINE  GRID 

COMMON  2T(30>,XI(60>,RO(10>,F(30,60,10).A(10>,B<10>,D(10),E<10), 

• EE(30),FF(30>,SV(30>,DPHI(10),BV{10),US(30,60>,UN(30,6C), 

• DNDS(30,10,2'.FAK(60>.FHW(60>.UTA(30>.UTB(30>,UTC(30>.BBV(30>, 

• BSA(30>.BSB(30).BSC(30  > , VSA( 30  > ,VSB( 30  > ,BAV{ 30  > ,LOP( 30  ) 

• .SONIC! 10) 

COMMON  DAF, DBF. DCF, AKLKE.RX, OX, CONST 

COMMON  H,DR0.DZT,TDZ.TPB.XIC,XIB,XID,TDR.AKLK,DR02,DZT2 
COMMON  0MCU,RHUB,RTIP,R2,RXE,RXH,UA,UB,TRDZ,ERA,ERB,CKLKE .SB, 

• BTR.FXB.FXI , EM2 , AMA . AM8 . DM . DKA . DKR , DO , DOS . TZ , TTZ , TZS , TRS , TOO 
COMMON  K.N.IDM,  KMX , LMX , NMX , KMXMl , KMXM2 . LMXMl , LMXM2 , NMXMl , 

• NMXM2,KTE0. IRX,  I TK . I TR , H RMX , KLEP , ITPR , KUP . KDN 
COMMON  ID(36) 

XIC  = 0.0 

XID  = OMCU 

XIB  = FXa*XID 

XII  = XIO  + FXI*XIO 

DO  = 2.0/FLOAT(KMX  + 1 > 

TOO  = 2.*DQ 
DOS  = DQ*DQ 
XIMP  = (XIB  + XI  I )/2. 

TALF  = ALOG( DQ/( 2.0-DQ) )/(  FXB*X I D-X IMP > 

ALF  = TALF/2. 

DO  1 K = 1 ,KMX 
Q = -1.0  + FLOAT! K >*DQ 
rAK!K)  = ALF*! 1 . - Q*Q > 

BKT  = DO*!  0.5<-FLOAT!  K ) >-l  . 

FHW!K)  = ALF*!1.  - 3KT*BKT) 

1 XI!K)  = XIMP  + ALOG! ! 1 .0  + Q)/! 1 ,0-0)  > / TALF 
4 OZT  = TPB/FLOAT! LMXMl ) 

DZT2  = DZT*DZT 
TDZ  = 2.0*0ZT 
DO  2 L = 1 .LMX 

2 ZT! L ) = FLOAT! L-1  )*DZT 

DRO  = ! RTIP-RHUB  )/FLOAT! NMXMl  ) 

TDR  = 2.*DRO 
DR02  = DRO*DRO 
DO  3 N = 1 ,NMX 

3 RO!N)  = FLOAT! N-1  )*DRO  + RHUB 
TZ  = DQ/DZT 

TZS  = TZ*TZ 
TTZ  = DQ./TDZ 
TRS  = D0S/DR02 

DAF  = !XI!KMX)-XI!KMXM1>)/!!XI!KMXM1>-XI!KMXM2>)*!XI!KMX) 

• -XI!KMXM2>>> 

DBF  = !XI!KMX>-XI!KMXM2>>/!!XI!KMXM1 >-XI!KMXM2>)*!XI!KMXMl > 

• -XI!KMX>)) 

DCF  = !2.*XI!KMX)-XI!KMXM1  )-XI!KMXM2>)/((XI!KMX>-XI!KMXM2>)* 

• ! XI! KMX  )-XI! KMXMl  ))  ) 

IF! IDM.EQ.2  ) GO  TO  13 
DO  12  N « 1 .NMX 

R2  = RO!N)*RO!N) 

E! N ) » 1.0  + R2 
D!N)  - !E!N)*E!N))/R2 
B!N)  = -E!N) 
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ISN 

0049 

12 

A<N)  = R2 

ISN 

0050 

GO  TO  14 

ISN 

0051 

13 

DO  15  N = l.NMX 

ISN 

0052 

R2  = RO(N)*RO(N> 

ISN 

0053 

E(N>  = 0.0 

ISN 

0054 

D(N>  = (1 .+R2)*( l.+R2>/R2 

ISN 

0055 

B<  N ) = -1  .-R2 

ISN 

0056 

A(  N ) = R2 

ISN 

0057 

15 

CONTINUE 

ISN 

0058 

14 

CONTINUE 

ISN 

0059 

DO  40  K = l.KMX 

ISN 

0060 

IF(XI(K).GT.XID)  GO  TO  41 

ISN 

0062 

40 

CONTINUE 

ISN 

0063 

41 

KTEO  = K 

ISN 

0064 

WRITE!  6,201  > XID 

ISN 

0065 

201 

FORMAT! /5X, 'THE  BLADES  LIE  BETWEEN  Z = 0 AND',F9.6) 

ISN 

0066 

WRITE!6,202) 

ISN 

0067 

202 

F0RMAT!//4X, 'K' ,7X, ' Z'.lOX,'  X/CA  '.SX.'RXE',  9X,'RXH',/) 

ISN 

0068 

DO  11  K = 1,KMX 

ISN 

0069 

PCT  = XI!K>/X1D 

ISN 

0070 

RX  = 1 . + ! RXE-1  . )*EXP! -!  ! XI ! K )-0. 5*XID >/XID )**2 ) 

ISN 

0071 

IF! ITPR.EQ.O  ) RX  = RXE 

ISN 

0073 

WRITE!  6,203  ) K , X I ! K ) , PCT , RX , RXH 

ISN 

0074 

1 1 

CONTINUE 

TSN 

0075 

203 

FORMAT! 3X,12,4X,F9.6,5X,F7.4,5X,F7.4,5X,F7.4> 

ISN 

0076 

WRITE!  6,204  ) 

ISN 

0077 

204 

FORMAT! //4X, 'L  ZETA',/) 

ISN 

0078 

DO  50  L = 1,LMX 

ISN 

0079 

WRITE!6,205)  L,ZT!L> 

ISN 

0080 

50 

CONTINUE 

ISN 

0081 

205 

FORMAT! 3X, I2,F9.4,8X,F9.6) 

ISN 

0082 

KHW  = INT! FLOAT! KMX )/2 . ) 

ISN 

0083 

WRITE!6.208)  KHW 

ISN 

0084 

208 

FORMAT! //4X, 'OPTIMUM  SHOCK  CAPTURING  OCCURS  WHEN  THE  GRID-SIZI 

• ' RATIO  RHO»!DELTA  ZETA)/!DELTA  Z)  » ! 1 +RHO*'*2  ) /RHO  ' , 

* /4X,'THE  TABLE  BELOW  LISTS  THESE  OPTIMUM  VALUES,  COMPARED 
'TO  THE  VALUES  ACTUALLY  USED  AT  K »',I2,) 

ISN 

0085 

WRITE!6,206)  KHW 

ISN 

0086 

206 

FORMAT! //4X, 'N  RHO ', 7X ,' US/WO  CRIT  ' , 5X,'  ARCTAN! RHO  ) , ' , 

* 6X,'M  REL '. 16X, 'GRID-SIZE  RATIO', 

/36X,  'DEGREES'  , 

* 26X, 'OPTIMUM  ACTUAL!AT  K = ',12,')',/) 

ISN 

0087 

DO  10  N = 1 ,NMX 

ISN 

0088 

BTR  = 1.0  - EM2*! 1 .O+RO! N )»RO! N ) ) 

ISN 

0089 

IF! SB.EQ.O.  ) GO  TO  51 

ISN 

0091 

USCRIT  = BTR/! SB«! 1 . +RO! N )*RO! N > ) > 

ISN 

0092 

51 

PSI  » 57.29578*ATAN! RO! N ) ) 

ISN 

0093 

SONIC!N)  = USCRIT 

ISN 

0094 

REL  = SORT!  EM2'*!  1 .+A!  N ) > ) 

ISN 

0095 

OPT  = ! 1 . +RO!  N )*RO! N ) )/PO! N ) 

ISN 

0096 

ACT  - RO!  N )’*0ZT»FHW!  KHW)/DQ 

ISN 

0097 

WRITE!  6,207  ) N,RO! N ) , USCR I T , PS  1 , REL ,OPT, ACT 

ISN 

0098 

10 

CONFINUb 

ISN 

0099 

207 

FORMAT!3X, I2,F9.4, IP£15.4,6X,0PF8.3,9X,F6.4,11X,FS.4,9X,F6.4> 

ISN 

0100 

RETURN 

ISN 

0101 

END 

B-8 


LEVEL  21.7  ( DEC  72  > 


OS/360  FORTRAN  H 


'ILER  OPTIONS  - NAME=  MAI N ,OPT=02 , L 1 NECNT=60 , S I2E=0000K , 

SOURCE, EBCDIC. NOLI  ST, NODECK, LOAD, MAP. NOEDI T, ID, XREF 
C 

SUBROUTINE  BVAL0( I > 

COMMON  ZT(30),XI(60).RO(10),F(30,60,10),A(10>,B(10),D(10>,E(10>, 

• EE( 30  ).FF<  30  >,SV( 30  >.DPHI( 10)  .BV( 10  ) ,US( 30,60>,UN{  30,60  >, 

• DNDS<  30 , 10,2  > ,FAK( 60  > ,FHW( 60 ) ,UTA( 30 ) ,UTB( 30  > ,UTC<  30 ) ,BBV( 30) , 

• BEA<.  30).BSB(30),BSC<30).VSA(30),VSB(30),BAV(30),LOP(30) 

• ,&ONIC(10) 

COMMON  DAF .DBF .DCF , AKLKE , RX , OX , CONST 

COMMON  H,DRO.OZT,TDZ,TP3,XIC,X1B,XID,TDR,AKLK,DR02,DZT2 
COMMON  0MCU.RHUB,RTIP,R2.RXE,RXH,IJA,UB,TRDZ,ERA,ERB,CKLKE,SB. 

• BTR.FXB.FXI , EM2 , AMA , AM3 , DM , DKA , DKR , DO , DOS , TZ .TTZ . TZS , TRS , TDQ 
COMMON  K.N.IDM,  KMX , LMX , NMX , KMXMl , KMXM2 , LMXMl , LMXM2 , NMXMl , 

• NMXM2,KTE0,  IRX,  I Tl<  , I TR  , I TRMX  , KLEP  , 1 TPR  . KUP  . KDN 
COMMON  ID(36) 

1F<  XI  ( K ) . LT. XIC  ) GO  TO  1 
IF(X!(K).LT.XID)  GO  TO  2 
GO  TO  3 

1 EE(  1 ) = 0. 

FF(  1 ) = FfLMX.K.N) 

GO  TO  10 

2 GO  TO  ( 100,200  ) , I 
200  EE<  1 ) = 0. 

KB  - K “ KLEP 

FFd)  = ."(LMX.K.N)  + DNDS(KB,N.l) 

GO  TO  10 

100  SLP  = DMOri K-KLEP ,N, 1 ) 

AE  = -BSA(  1 )*2.*FAK<K)/RX  + 2 . ‘BSBI 1 )*FHW( K-1 ) -3.5*AKLKE 

• - BSC(  1 )»FHW( K-1  ) +4.*DKA 
BE  - -8.*0KA  + 4.*AKLKE 

CE  = BSA( 1 )*( FHW( K )*F( 1 ,K+1 ,N )-2 . *FAK( K )*OX*UTA( 1 ) + FHW( K-1 ) 

• ‘Fd.K-l.N))  + BSB(  1 )»(-FHW(K-l  )*(UTA(  1 ) + F(  1,K-1.N)) 

• -FHW( K-2 )•( F( 1 ,K-1 .N  )-UTC< 1 ) ) ) + BSC ( 1 ) * ( FHW( K )•( F ( 1 , K+ 1 , N > 

• -UTA(  1 ) ) + FHW< K-1  )*F( 1 ,K-1  ,N  ) ) 

• -DKA*(UTC(3)-UTA(3)-8.»UTB(2)  + 8.*UTB( 1 )-UTC( 1 )-3.*UTA( 1 )) 

• +AKLKE»(-.5*UTA(3)-3.*RO(N)*DZT*SLP-A<N)*DZT*UB»(UTC< 1 ) 

• -6.*UTB( 1 ) + 3.«UTA( 1 ) + 2.*F( 1 , K+ 1 , N ) )*FAK( K ) ) 

• +DI;R*(ERA*F(1,K,N+1)+ERB*F(1,K.N-1)-2.»F(1,K.N)) 

EE(  1 ) = -BE/AE 

12  FF(  1 ) = -CE/AE 
GO  TO  10 

3 EEd  ) = 0.0 

Z = F( LMX.K.N) 

6 FF<  1 ) = Z ♦ DPHKN  ) 

10  RETURN 
END 


LEVEL 


ISN 

ISN 


ISM 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 


ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 


ISN 


21.7  < DEC  72  > 


OS/360  FORTRAN  H 


COMPILER  OPTIONS 
C 


NAME=  MAIN,OPT=02,LINECNT=60,SIZE=0000K, 

SOURCE, EBCDIC, NOLI  ST, NODECK, LOAD, MAP, NOED IT, ID, XREF 


0002 

0003 


0004 

0005 

0006 

0007 

0008 

0009 

0010 
001  1 
0012 
0013 
0015 

0017 

0018 

0019 

0020 
0022 

0023 

0024 

0025 


0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 


0036 


SUBROUTINE  BVAL 1 ( I > 

COMMON  ZT( 30 ),XI<60  >,RO< 10),F(  30,60, 10),A( 10>,B( 10),0( 10>,E< 10>, 

* EE(30),FF<30),SV(30),DPH1( 10),BV{ 10),US(30,60>,UN<30,60>, 

* DNOS(30,10,2),FAK<60>,FHW(60>,UTA(30),UTB(30>,UTC<30>,BBV{30>, 

* BSA( 30  > ,BSB( 30 ) , BSC ( 30  > , VSA<  30 ) , VSB( 30 ) ,BAV{ 30  > ,LOP( 30  ) 

* ,S0NIC(10) 

COMMON  DAF ,DBF ,DCF , AKLKE , RX , OX , CONST 

COMMON  H,DRO,OZT, TDZ,TPB,XIC,XIB,XID,TDR,AKLK,DR02,DZT2 
COMMON  OMCU.RHUe ,RTIP.R2,RXE,RXH,UA,UB,TRD2,ERA,ERB,CKLKE ,SB, 

* PTR , F y B , F X 1 , EM2 , AMA , AMB , DM . OKA , OKR , DQ , DOS . TZ , TTZ ,TZS ,TRS , TDQ 
COMMON  K,N,IDM,  KMX , LMX , NMX . KMXMl , KMXM2 , LMXMl , LMXM2 , NMXMl , 

* NMXM2,KTE0, IRX,  ITK , I TR , ITRMX , KL E P , I TPR , KUP , KDN 
COMMON  id; 36) 

FPA  = F( 2,K-1 ,N  > 

FPB  = F( 2,  K ,N  > 

FPC  = F(2,K+1 ,N) 

DKL  = 0. 

IF(XI<K).LT.XIC)  GO  TO  1 
IF(yi(K>.LT.XID>  GO  TO  2 
FPA  = FPA  - DPHHN) 

FPB  = FPB  - DPHKN) 

FPC  = FPC  - DPHKN) 

IF( I0M.EQ.2  ) GO  70  1 
DKL  = ( 1 . +RO<  N )»RO( N ) )• 

* 0.5*TRS*{ERA*DPHI{N  + 1 ) + ERB»DPHI ( N-1  ) -2. ‘DPHl ( N ) >/FAK( K ) 

1 L = LMX 

BKL  = -2. ‘AKLKE  - 2 . *FAK( K )»BSA( L )/RX  + 2 . *FHW(  K- 1 ) ‘BSB;  L ) 

* - BSC( L )*FHW( K-1  ) 

DKL  = DKL  - BSA(  L )•( FHWl K )»F< L ,K+1 ,N )-2.*FAK< K>*OX*UTA{ L ) 

* + FHW(K-1 )*F(L,K-1,N))  + BSB ( L ) •< FHW< K- 1 >*(  UTA< L ) + F ( L , K- 1 , N ) ) 

* + FHW;  K-2 )•(  F( L ,K-1 ,N  )-UTC( L ) ) ) - BSC(  L )*< FHW( K )*( F< L . K + 1 , N ) 

* -UTA(L))  + FHW( K-1  )»F( L ,K-1 ,N  ) ) 

* +2  . ‘DKA*;  -FPA  + F(  I.  ,K-1  ,N)  + FPB-2.*'F(L,K,N)  + F(L-1,K,N) 

* +F(L,K+1,N)-F(L-1.K+1,N)) 

* - DKR*< ERA*F( L ,K,N+1  ) + ERB*F( L ,K,N-1  )-2.*UTA( L ) ) 

61  AKL  = AKLKE 

CKL  = CKLKE 

25  SV(LMX)  = (DKL-AKL*FPB  -CKL*FF ( LMXMl  ))/( BKL +CKL«EE ( LMXMl  ) > 

GO  TO  10 

2 CONTINUE 

GO  TO  ( 100,200), I 
200  SLP  ^ DNDS<K-KLEP,N,2) 

AE  = -2.*FAK(K)*BSAILMX)/RX  + 2 . ‘BSBl LMX  )*F HVM  K- 1 > - 7.*AKLKE/2. 

* -BSC( LMX )*FHW( K-1  ) + DKA 
BE  = -8.*DKA  + 4. ‘AKLKE 

CE  = BSA(  LMX )*< FHW( K )‘F( LMX,K+1 ,N )-2 .‘FAK; K )*UTA( LMX  )‘OX 
‘ + FHW; K-l  )‘F( LMX ,K-1 ,N  ) ) + BSB< LMX )‘( -FHW< K-1 )*( UTA< LMX > 

* + F( LMX,K-1 ,N  ) ) - FHW(K-2)‘(F(LMX,K-1,N)-UTC<LMX))) 

‘ + BSC<  LMX )‘(  FHW;  K)‘(  F;LMX.K+1 ,N )-UTA( LMX ) ) + FHW{K-l ) ‘F < LMX , K- 1 , N ) ) 

‘ -DKA‘(F(LMXM2,K  + 2,N )-UTA(LMXM2  )-8.‘F( LMXMl , K+1 ,N)  + 4.‘( 

‘ F(  LMX,i<+l  ,N  )-UTB(  LMX  ) ) + UTC<LMX)) 

‘ +AKLKE‘(-F(LMXM2,K,N)-6.‘DZT‘SLP  - 1 1 . ‘UTA< 1 ) + 1 8 . * 

‘ UTA( 2 )-9. ‘UTA( 3 ) + 2.‘UTA( 4 ) )/2. 

‘ +DKR‘(ERA‘F(LMX,K,N  + 1)  + ERB*F(LMX,K,N-1 ) -2 . ‘F < LMX , K , N ) ) 

GO  TO  11 


B-10 


ISN  0037 
ISN  0038 
ISN  0039 


100  SLP  = DNDS<K-KLEP,N,2) 

L = LMX 

AE  = -2 . *FAK( K >*BSA( LMX  )/RX  + 2 . ‘ESB ( LMX >*FHW( K- 1 > 

• -BSC( LMX )*FHW< K-1 ) + OKA  - 3.5*AKLKE 

ISN  0040  BE  = 4.*A'<LKE  -8."D|'A 

ISN  0041  CE  = 3SA<LMX>''(FHV/(K)*F(LMX.K+1  ,N)-2.*FAK(K)*UTA(LMX)»0X 

• +FHW<  K-l )*F<  LMX,K-1 ,M  > > + BSB( LMX )*( -FHW( K-1 )•( UTAC  LMX ) + F( LMX, 

• K-1  ,N  ) )-FHV/<  K-2  )•(  F(  LMX.K-1  ,N  )-UTC(  LMX  > ) ) + BSC(  LMX  >•(  FHW(  K >* 

• < F<  LMX .K+1 .N )-UTA( LMX  >>  + FHW<K-l)*F<LMX.K-l,N)> 

• -DKA*(F<  LMXM2,K  + 2.N  >-UTA(LMXM2>-0.*F(LMXMl ,K+1 ,N)  + 4.»F<LMX, 

• K+1  ,N  >-4  . ‘UTSf  LMX  > + UTC(  LMX ) )+AKLKE«( -0.5*UTA<  LMXM2 ) + 3 . ‘ROl N ) 

• *DZT*SLP+A(  N > *021*118 •{  UTC<  LMX  )-6  . *UTB ( LMX  ) + 3 . *UTA(  LMX  > + 2.* 

• F( LMX .K+1 ,N  > )*FAK( K ) ) 

• +DKR*(  ERA*F(  LMX,K.N  + 1 ) + E RB*F  ( LMX , K . N- 1 )-2.*F(LMX,K,N)> 

ISN  0042  1 1 SV(LMX>  = -< CE + BE *r F ( LMXM 1 ) > / ' AE + 3E *EE ( LMXMl  ) ) 

ISN  0043  10  RETURN 

ISN  0044  END 


LEVEL  21.7  < DEC  72  > 


OS/360  FORTRAN  H 


ISN 

ISN 


ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISM 

ISN 


COMPILER  OPTIONS  - NAME=  MAI N ,OPT=02 . L I NECNT=60 , S I ZE=0000K. 

SOURCE.EBCDIC.NOLIST.NODECK.LOAD.MAP.NOEOIT, ID.XREF 
C 


0002 

0003 


0004 

0005 

0006 

0007 

0008 

0009 

0010 
0012 
0014 
0016 
0018 


SUBROUTINE  OUTPUT! lOP , I BC  ) 

COMMON  ZT(30>,XI(60>,RO(IO),F(30,60,10>.A<10),B(10),D(10>,E<10>, 

• EE(30),FF(30),SV(30),DPHI(10>,BV(10),US(30,60>,UN(30.60), 

• DNDS<30,10,2 ) ,FAK( 60 ) .FHW( 60 ) ,UTA( 30 ) .UTB( 30  > ,UTC<  30  ) ,BBV<  30  ) , 

• BSA(30),BSB(30),BSC(30 ) ,VSA( 30 ) ,VSB( 30  ) ,BAV( 30  ) .LOP( 30  > 

• ,SONIC<10> 

COMMON  DAF .DBF ,DCF , AKL KE , RX . OX , CONST 

COMMON  H . DRO , DZT , TDZ ,TPB,XIC,XIB,XID,TDR, AKLK , DR02 . DZT2 
COMMON  OMCU ,RHUB ,RTIP , R2 . RXE , RXH , UA . UB .TROZ . ERA , ERB .CKLKE , SB , 

• BTR.FXB.FXl , EM2 , AMA , AMB , DM . DKA . DKR , OQ , DOS , TZ , TTZ , TZS , TRS , TDQ 
COMMON  K.N.IDM,  KMX . LMX , NMX , KMXMl , KMXM2 , LMXMl . LMXM2 , NMXMl , 

• NMXM2,KTE0, IRX,  I TK , I TR , I TRMX , KL EP . I TPR , KUP  . KDN 

COMMON  ID(36> 

DIMENSION  UR(30,60> 

IF< IOP.EQ.2)  GO  TO  2 
IF( lOP .EQ. 3 ) GO  TO  20 
IF( IOP.EQ.4)  GO  TO  20 
IF< lOP .EQ.5  > GO  TO  90 
IF( IOP.EQ.6  ) GO  TO  96 


C 

C lOP  = 1 GIVES  VALUES  OF  THE  POTENTIAL  AMD  ALL  VELOCITY  COMPONENTS 
C AT  all  grid  POINTS. 

C lOP  2 GIVES  ONLY  VELOCITY  COMPONENTS,  AND  ONLY  AT  ZETA  = 0.  AND  2.*PI/B 

C lOP  = 3 IS  THE  SAME  AS  lOP  = 1 , BUT  WITH  VALUES  OF  THE  POTENTIAL 

C OMITTED. 

C lOP  = 4 GIVES  MACH  NUMBER  CONTOURS. 

C lOP  = 5 CALCULATES  PERFORMANCE  DATA 

C lOP  = 6 GIVES  THE  LOCAL  MACH  NUMBER  AT  ALL  POINTS 

C lOP  = 7 PRINTS  ONLY  THE  POTENTIAL  (AT  ALL  POINTS) 

C 

C PRINT  PHI 
C 


ISN 

0020 

ISN 

0021 

ISN 

0022 

ISN 

0023 

250 

ISN 

0024 

10 

ISN 

0025 

200 

ISN 

0026 

ISN 

0027 

ISN 

0028 

80 

ISN 

0029 

251 

ISN 

0030 

ISN 

0031 

ISN 

0033 

ISN 

0034 

ISN 

0036 

ISN 

0037 

20 

ISN 

0938 

ITR, ITK,N,RO(  N ) 

OF  THE  POTENTIAL  AFTER  ITR 


LA  ’ 1 
LB  = 10 
WRITE!  6,250  > 

FORMAT! IHl , 

* //3X, 'VALUES 

‘ ' . ITK  = ■ ,14. 

• ' AT  RHO!  ',12,')  = ',F7.4) 
WRITE (6,200  ) (LOPCl  ),L  = LA,LB) 
FORMAT! /4X.  'K  L= ' , I 3 , 9 I 1 2 , // ) 

DO  30  K = l.KMX 

WRITE (6, 251)  K.(F!L,K.H),L  = LA,LB) 

CONTIHUt 

FORMAT! 15, 1P10E12.3) 

LA  = LA  + 10 

IF! LA. GT  LMX  ) GO  TO  20 

LB  = LB  + 10 

IF! LB  .GT.LMX ) LB  = LMX 

GO  TO  10 

CONTINUE 

IF! lOP .EQ.7  ) RETURN 


M4. 


C 

C CALCULATE  US  AND  UN 
C 


B-12 


ISN  0040 
ISN  0041 
ISN  0042 
ISN  0043 

ISN  0044 
ISN  0045 
ISN  0046 
ISN  0047 
ISN  0048 
ISN  0050 
ISN  0052 

ISN  0053 

ISN  0054 
ISN  0055 
ISN  0056 
ISN  0057 
ISN  0058 


00  60  K = 2,I<MXM1 
DO  6!  L = 2,LM;!M1 

US(L,K!  UB*,-AK<K)»(F(L,K+1,N>  - F(L,K-1.N)> 

61  UN(L.K)  = ( ( F( L + 1 ,K-1 ,N  )-F( L ,K-1 ,N > ) + ( F( L ,K+1 ,N  )-F<L-l ,K+1 ,N > > ) 
• /T1D2  - UA*US(L,K> 

60  CONTINUE 

DO  63  ■<  = 2,KMXM1 

USd.K)  = IJB*FAK(K)»(F(1,K  + 1.N)-F(1.|!:-1.N>) 

USd-HX.Kl  = U8*FAK( K )*f F( LMX.K+1 .N >-F< LMX,K-1 ,N  ) ) 

IFIK.LE.KLEP  ) CD  TO  62 
IF< K.CE .‘TEO > GO  TO  62 

UN<1,K)  - ( -3.  *F(  1 ,K,N  ) + 4 .'•■•FI  2,K-1  ,N  >-F(  3.K-2.N  ) 


* +F( 1 ,K+1 ,N )-F( 1 ,K-I ,N > )/TRDZ  - UA»U5(1,K) 

UN(LMX,K)  = ( 3 .*F( LMX ,K,N >-4 .*F( LMXMI .K+1 ,N >+F< LMXM2,K+2,N > 

* + F(LMX,K+1 ,N)-F(LMX,K-1 ,N)  )/TRDZ  - UA*US(LMX,K> 

GO  TO  63 

62  FZTA  = ( F<  2 ,K-1  ,N  )--F(  1 ,K-1  ,N  ) + F(  LMX  ,K+I  ,N  )-F<  LMXMl  ,K+1  ,N  > >/TRDZ 
UNI  I ,K  > = FZTA  - UA*US( 1 , K ) 

UNI LMX,  K>  = FZTA  - UA*USILMX,K> 

63  CONTINUE 


C 

C FIND  MACH  NUMBER  CONTOURS 
C 


ISN 

0059 

IFI iCP . NE . 4 > GO  TO  11 

ISN 

0061 

TGTM  = AMA 

ISN 

0062 

RELSO  = EM2»( 1 .+AIN)  > 

ISN 

0063 

GMAP  ’ SE/EM2 

ISN 

0064 

REL  = SQRriRELSQ) 

ISN 

0065 

CF  = GMAP/2. 

ISN 

0066 

504  CON  1 INUE 

ISN 

0067 

CALL  CLEAR  I URI 1 . 1 ) , URI  30 . 60  ) > 

ISN 

0063 

CALL  CLEAR  I UNI  1 , 1 ) , UNI  30 , 60  ) ) 

C 

C THE  FOLLOWING  FORMULA  WAS  USED  IN  CALCULATING  THE  MACH  NUMBER 
C CONTOURS  IN  AIAA  PAPER  77-199 
C TGT  = I ITGTM»TGTM/REl.SQ)-l  . >/GMAP 

C 

C THE  FOLLOWING  IS  THE  LINEARIZED  FORMULA 
C 


ISN 

0069 

TGT  =-I 1 .-TGTM/REL )/CF 

ISN 

0070 

KWMX  = 0 

ISN 

0071 

DO  300  1.  = l.LMX 

ISN 

0072 

SGN  = lUSIL.KUP)  - TGT)/ABSIUSIL, 

ISN 

0073 

KW  = 1 

ISN 

0074 

DO  501  K = KUP.KDN 

ISN 

0075 

IFI SGN.GT.O.  ) GO  TO  502 

ISN 

0077 

IFIIUSIL.KI-TGTI.GT.O. ) GO  TO  503 

ISN 

0079 

GO  TO  501 

ISN 

0080 

502 

IFIIUSIL.IO-TGTI.LT.O.  ) GO  TO  503 

ISN 

0082 

GO  TO  501 

ISN 

0033 

503 

X2  = XIIK) 

ISN 

0084 

XI  = XII K-1  ) 

ISN 

0085 

VZ  = U3I L ,K> 

ISN 

0086 

Y1  = U3IL,K-1 > 

ISN 

0087 

Z = XI  + I X2-X1  )*(TGT-Y1  )/( Y2-VI > 

ISN 

0083 

URIL.KW)  =■  Z/XTD 

ISN 

0089 

UNIL.KW)  = ROIN)*IZTIL)+Z>/OMCU 

B-13 


ISN  0090 
ISN  0092 
ISN  0093 
ISN  0094 
ISN  0095 
ISN  0096 
ISN  0093 
ISN  0099 


ISN  0100 
ISN  0101 
ISN  0102 
ISN  0104 
ISN  0105 
ISN  0106 
ISN  0107 
ISN  0108 
ISN  0109 
ISN  0110 
ISN  0111 
ISN  0112 
ISN  0113 
ISN  0114 
ISN  0115 
ISN  0116 
ISN  0113 
ISN  0119 
ISN  0121 
ISN  0122 
ISN  0123 
ISN  0125 


ISN  0126 
ISN  0127 
ISN  0128 
ISN  0129 

ISN  0130 
ISN  0131 
ISN  0132 
ISN  0133 
ISN  0134 
ISM  0135 
ISN  0137 
ISN  0138 
ISN  0140 


ISN  0141 
ISN  0142 
ISN  0143 
iSN  0144 

ISN  0145 


IF<  KW.GT.KWMX  ) KWMX  = KVI 
KW  = KW  + 1 
SGN  = -SGN 
501  CONTINUE 
500  CONTINUE 

IF( KWMX . EQ. 0 > GO  TO  510 
WRITE(6,205)  N,RO(N).TGTM 

205  FORMAT! /5X. 'FOR  RO('.!2.’l  = ’,F7,4,'  THE  LOCAL  MACH  NUMBER' . 

* ' EQUALS'  ,F7. 4,  ' AT  THE  FOLLOWING  VALUES  OF  X/CA  AND  ', 

* ' R*THETA/CA ' , / ) 

LA  = 1 

L6  = 20 

IF(LMX.LT.20  ) LB  = LMX 

511  WRITE(  6.230  ) ( LOP ( L ) . L = LA , LB  ) 

230  FORMAT! / 1 X KW  L = ' , I 2 , 1 9 1 6 , // > 

231  FORMAT!  rX,I2,20F6.2) 

DO  512  Kl/  = 1 .KWMX 

WRITE! 6 , 231  ) KW, ! UR! L ,KW) . L = LA,LB  > 

WRITE!  6, 206  > ! U N ! L . KW ) , L = L A , L B ) 

WRITE! 6 ,207  > 

512  CONTINUE 
520  CONTINUE 

206  FORMAT!  31<,20F6.2  > 

207  FORMAT! IX) 

LA  = LA  + 10 

IF! LA.GT.LMX  > GO  TO  510 
LB  = LB  ••  10 
IF!LB.GT.LMX)  LB^LMX 
GO  TO  511 

510  TGTM  = TGFM  + DM 

IF! TGTM.LE .AMB  ) GO  TO  504 
RETURN 
C 

C PRINT  US 
C 

11  LA  - 1 
LB  = 10 

WRITE!  6,2  53  ) I TR . ITK , N . RO! N ) .SONIC! N ) 

253  FORMAT! ///3X, 'VALUES  OF  US/WO  AFTER  ITR  - ',14,',  ITK  = ',14, 

• ' AT  RHO!'.I2,')  = ' , F7. 4 , lOX, ' US/WO  CRIT  = ',1PE12.3) 

21  WRITE!  6,200>  ! LOP!  I ) , I =LA  , LB  ) 

DO  31  K = 2,KM)!M1 

WRITE!6,251>  K,!US!L,K).L=LA,LB) 

81  CONTINUE 

LA  = LA  + 10 
IF! LA.GT.LMX)  GO  TO  40 
LB  = LB  + 10 
IF! LB.GT. LMX > LB  = LMX 
GO  10  21 
C 

C PRINT  UN 
C 

40  WR1TE!  6,254  > I TR , 1 TK , N , RO! N ) 

LA  - 1 

LB  = 10 

254  FORMATI///3X, 'VALUES  OF  UN/WO  AFTER  ITR  » ',14,',  ITK  - *,I4, 

• ' AT  RHOI '.12,')  = ‘,F7.4> 

41  WRITE!  6.200  > ! LOP ! 1 ) . I = LA , L B ) 
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ISN  0146 
ISN  0147 
ISN  0143 
ISN  0149 
ISN  0150 
ISN  0152 
ISN  0153 
ISN  0155 
ISN  0156 
ISN  0158 


ISN  0174 
ISN  0175 


ISN  0189 
ISN  0190 


ISN  0191 
ISN  0192 
ISN  0193 
ISN  0194 
ISN  0196 


DO  82  K = Z.KMXMl 

WRITE! 6,251  ) K , < UN( L . K ) , L = LA, LB > 

82  CONTINUE 

LA  = LA  + 10 
IFILA.GT.LMX)  GO  TO  31 
LB  = LB  + 10 
IF( LB.GT.LMX ) LB  = LMX 
GO  TO  41 

31  IF< IDM.EQ. 2 ) WRITE! 6,210 > 

210  FORMAT! //lOX, 'NOTE:  THE  FOLLOWING  VALUES  OF  V RADIAL  ARE  BASED', 
• ' ON  THE  2D  STRIP-THEORY  APPROXIMATION') 

CALCULATE  V RADIAL 


ISN 

0159 

IF!  N. 

EQ. 

1 ) 

GO  TO  52 

ISN 

0161 

IF!  N. 

EQ. 

NMX  ) 

GO  TO  52 

ISN 

0163 

DO  50 

K 

s 

2, 

KMXMl 

ISN 

0164 

DO  51 

L 

= 

1 , 

LMX 

ISN 

0165 

UR!  L , 

K ) 

s 

! F 

< L , K , N + 1 ) 

ISN 

0166 

51 

CONTI 

NUE 

ISN 

0167 

50 

CONTI 

NUE 

ISN 

0168 

GO  TO 

53 

ISN 

0169 

52 

DO  64 

K 

s 

2, 

KMXMl 

ISN 

C170 

DO  55 

L 

= 

1 , 

LMX 

ISN 

0171 

UR!  L , 

K) 

s 

0 . 

ISN 

0172 

55 

CONTI 

NUE 

ISN 

0173 

54 

CONTI 

NUE 

PRINT  V RADIAL 

53  WRITE!  6,255  > ITR , ITK , N , RO! N ) 

255  FORMAT! ///3X, 'VALUES  OF  V RADIAL/U  INFINITY', 

-*  'AFTER  ITR  = ',14,',  ITK  = ',14, 

* ' AT  RHO!  ' ,12,  ' ) = ' ,F7.4) 

LA  = 1 

LB  = 10 

50  WRITE!  6,200  ) ! LOP ! I ) , ! =LA , L B ) 

DO  83  K = 2,I<MXM1 

WRITE! 6,251  ) K , ! UR! L , K ) . L = LA , LB ) 

83  CONTINUE 

LA  = LA  + 10 
IF!  l.A.GT.LMX)  GO  TO  30 
LB  = LB  + 10 
IFILB.GT.LMX)  LB  = LMX 
GO  TO  56 

CALCULATE  AND  PRINT  VELOCITIES  ON  ZETA  = 0 AND  ZETAl 

2 WRITE!6,201)  I TR , ITK , N , RO! N ) , SON IC! N ) 

201  FORMAT!//4X, 'SURFACE  VELOCITIES  AFTER  ITR  • ',14,',  ITK  - ',14, 

• ' AT  RHO!', 12,')  = ' ,F7.4, lOX, 'US/WO  CRIT  - ',1PE12.3, 

* /4X,  ' K'  ,4X,  'US! 1 ,K  ) ' ,4X, 'US!LMX,K)  ' ,4X, 'UN! 1 ,K)' ,4X, 

• 'UN! LMX,K)  ' ,4X, 'UR! 1 ,K) ' ,4X, 'UR! LMX,K) ' ,/ ) 

DO  65  K = 2,KMXM1 

US!1,K)  » UB*FAK!K)*!F! 1 ,K+1 ,N)-F! 1 ,K-1 ,N) ) 

US!LMX,K)  • U3*FAK!K)«!F!LMX,K+1 ,N)-F!LMX,K-1 ,N) ) 

IF!K.LE.KLEP)  GO  TO  64 
IF!K.GE.KTEO)  GO  TO  64 
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ISN 

0198 

ISN 

0199 

ISN 

0200 

ISN 

0201 

ISN 

0202 

ISN 

0203 

ISN 

0204 

ISN 

0206 

ISN 

0208 

ISN 

0209 

ISN 

0210 

ISN 

0211 

ISN 

0212 

ISN 

0213 

ISN 

0214 

ISN 

0215 

ISN 

0216 

ISN 

0217 

ISN 

0218 

ISN 

0219 

ISN 

0220 

ISN 

0221 

ISN 

0222 

ISN 

0223 

ISN 

0224 

ISN 

0225 

ISN 

0226 

ISN 

0227 

ISN 

0228 

ISN 

0229 

ISN 

0230 

ISN 

0231 

ISN 

0232 

ISN 

u233 

ISN 

0234 

ISN 

0235 

ISN 

0236 

UNd.K)  = ( -3.*F(  1 ,K,N  ) + 4.*F<  2.K-1  ,N  )-F(  3.K-2,N) 

* +F( 1 ,K+1 ,N )-F{ 1 ,K-1 ,N  ) >/TRDZ  - UA*US{1,K) 

UN<LMX,K)  = ( 3. *F( LMX.K.N >-4 . *F< LMXMl ,K+1 ,N  ) + F( LMXM2,K  + 2,N  ) 

* + F(LMX.K+1,N>-F(LMX.K-1,N))/TRDZ  - UA*US<LMX,K) 

GO  TO  66 

64  FZTA  = ( F(  2,K-1 , N >-F( 1 .K-1 ,N ) + F( LMX ,K+1 ,N >-F( LMXMl ,K+1 ,N > >/TROZ 
UN( 1 , K ) = FZTA  - UA*US( 1 ,K> 

UN(LMX,K>  = FZTA  - UA’*US(  LMX  , K ) 

66  IF<N.EQ.  1 > GO  TO  70 
1F( N.EQ. NMX  > GO  TO  70 

URd.K)  = (Fd.K.N+l>  - Fd  ,K,N-1  > )/TDR 
UR(LMX,K)  = ( F( LMX.K, N+1  ) - F ( LMX , K , N- 1 ) ) /TDR 
GO  TO  72 

70  UR( 1 ,K>  = 0.0 
UR(LMX,K>  = 0.0 

72  WR1TE(6,251)  K , US(  1 , K ) , US ( LMX , K ) , UN (1 , K > , UN( LMX , K ) , 

* UR( 1 ,K ) , UR( LMX ,K ) 

65  CONTINUE 
30  CONTINUE 

K1  = KTEO 
K2  = KTEO  - 1 
K3  = KTEO  - 2 
XI  = XI( K1  > 

X2  = XI  ( K2  > 

X3  = XI(K3) 

Y!  = US(  1 ,K1  ) 

V2  = US(  1 ,K2  > 

V3  = US( 1 ,K3  ) 

X = XID 

Z = ( X-X2  )*( X-X3  )*Yl/( ( X1-X2  )*( X1-X3  ) > 

* +(X-X1>*{X-X3>*V2/<(X2-X1)*(X2-X3)) 

* +(X-X1)»<X-X2>»Y3/((X3-X1)*(X3-X2)) 

U1  = Z 

VI  = US(LMX,K1 > 

Y2  = US(LMX,K2> 

Y3  = US(LMX,K3> 

Z = ( X-X2 >*( X-X3  )*Yl/( < X1-X2 >*( X1-X3 > > 

* +(X-';i)»(X-X3)*Y2/{(X2-Xl)*(X2-X3)) 

* •KX-Xi)*(X-X2>*Y3/((X3-Xl>*(X3-X2)) 

U2  = Z 

DU  = U1  - U2 

WRITE(  6.202  ) U 1 , U2  , DU , N , DPH I < N > 

202  FORMAT!  /5X , ' TRA I L I NG-E DGE  VELOCITIES  ARE  US( 1 ,TE ) = '.IPEIO.S, 

* ' US(LMX,TE)  = ',E10.3,'  DELTA  US  = ',E10.3, 

* ' DPHK  ',12.')  = '.E10.3) 

IF ( IBC  .NE .2  ) RETURN 

C 

C CALCULATE  THE  BLADE  SHAPE 
C 


ISN 

0238 

IF<  I0P.NE.6) 

GO 

ISN 

0240 

K1 

= 

KTEO 

ISN 

0241 

K2 

= 

KTEO  - 

1 

ISN 

0242 

K3 

a 

KTEO  - 

2 

ISN 

0243 

XI 

= 

Xl( K1 > 

ISN 

0244 

X2 

a 

Xi  ( K2  > 

ISN 

0245 

X3 

a 

XI(K3  ) 

ISN 

0246 

KO 

a 

KLEP  - 

3 

ISN 

0247 

K4 

a 

KTEO  + 

3 
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ISN  0248 
ISN  0249 
ISN  0250 
ISN  0251 
ISN  0253 
ISN  0255 

ISN  0256 

ISN  0257 
ISN  0258 
ISN  0259 
ISN  0260 
ISN  0261 
ISN  0262 
ISN  0263 
ISN  0264 
ISN  0265 
ISN  0266 
ISN  0267 


ISN  0268 
ISN  0269 
ISN  0270 
ISN  0271 
ISN  0272 


ISN  0273 
ISN  0274 
ISN  0275 
ISN  0276 
ISN  0277 
ISN  0278 
ISN  0279 
ISN  0280 
ISN  0231 
ISN  0282 
ISN  0233 
ISN  0284 


ISN  0285 
ISN  0286 
ISN  0287 
ISN  0288 
ISN  0289 


ISN  0290 
ISN  0291 
ISN  0292 
ISN  0293 
ISN  0294 
ISN  0295 
ISN  0296 
ISN  0297 


DO  23  K = K0,(<:4 

US<1,K)  = UB*FAK(<)*(F(1,K+1,N)-F(1,K-1,N>) 

USILMX.K)  = UB*FAK(K>»(F(LMX,K+! ,N>-F(LMX,K-1,N>) 

IFIK.LE.KLEP)  GO  TO  24 
IF(K.GE.KTEO)  GO  TO  24 

UN(1,K)  = (-3.*F(1.K,N>+4.*F(2,K-1.N)-F<3,K-2,N) 

• +F< 1 ,K+1 ,N)-F( 1 ,K-1 ,N)  )/TRDZ  - UA*US<1,K> 

UNILMX.K)  = ( 3.*F( LMX.K.N  1-4 .‘FI LMXMl .K+1 ,N > + F( LMXM2,K  + 2,N > 

• + F( LMX ,K+1 ,N )-F( LMX,K-1 ,N  ) >/TRDZ  - UA*US(LMX,K) 

GO  TO  23 

24  FZTA  = < F< 2.K-1 ,N )-F< 1 ,N  ) + F< LMX,K+1 ,N  )-F( LMXMl ,K+1 ,N > )/TRDZ 
UN(I.K)  = FZTA  - UA*US<1,I<) 

UNILMX.K)  = FZTA  - UA^USI LMX , K ) 

23  CONTINUE 

22  ORT  = SQRT< 1 . +R0( N >*R0( N ) )/OMCU 
Y1  = UNI  1 .K1  > 

Y2  = UNI  1 ,K2  ) 

Y3  = UNI  1 .K3  ) 

X = XIO 

Z = IX-X2)»IX-X3>*Y1/I IX1-X2)*IX1-X3)) 

• +IX-X1)»IX-X3)»Y2/IIX2-X1)*I X2-X3  ) ) 

• +I X-Xl  )•! X-X2  >»Y3/I I X3-X1  )*IX3-X2  > > 

UTEU  = Z 

Y1  = UNI LMX.Kl  ) 

YZ  = UNILMX,K2) 

Y3  = UNILMX.K3) 

Z = I ;;-X2  )•!  X-X3  >*YI/I  I X1-X2  >*I  X1-X3  ) ) 

• +IX-X1)*IX-X3>»Y2/IIX2-X1)*IX2-X3)) 

• +IX-X1)*IX-X2)»Y3/II X3-X1 >*I X3-X2  ) ) 

UTEL  = Z 

K1  = KEEP  + 1 
K2  = KEEP  + 2 
K3  = KEEP  + 3 
Y1  = UNI  1 ,K1  ) 

Y2  •■=  UNI  1 ,K2) 

Y3  = UNI  1 ,K3  > 

XI  = XIIKl  > 

X2  = XIIK2> 

X3  = XI I K3  ) 

X = 0. 

Z - IX-X2)*IX-X3)*Y1/IIX1-X2)»IX1-X3)) 

• H X-Xl >•! X-X3 )»Y2/I I X2-X1  I”! X2-X3  > ) 

• + 1 X-Xl >*I X-X2)»Y3/I IX3-X1  l-IXS-XZ  ) ) 

UEEU  = Z 

Y1  = UNIEMX.Kl  ) 

Y2  = UNIEMX,K2) 

Y3  = UNIEMX,K3) 

Z = IX-X2)*IX-X3>*Y1/IIX1-X2>»IX1-X3)) 

• + 1 X-Xl  )•!  X-X3  )*Y2/I  IX2-Xn»IX2-X3>> 

• ♦IX-Xl )•! X-X2 )*Y3/I I X3-X1 >»IX3-X2)> 

UEEE  = Z 

USIl.KI)  > 0.5*XIIK1 >*IUNI l.Kl >+UEEU) 

USIEMX.Kl)  = 0.5*XlIKI)*IUHIEMX.Kn+ULEL) 

KA  = KEFP+2 
KB  = KTEO-1 
DO  75  K •=  KA.KB 

USIl.K)  = USIl,K-n  + 0.5»IXIIK)-XlIK-l))*«UN(l.K>+UNU,K-l>> 
USIEMX.K)  - USIEMX.K-l  ) + 0.5*IXIIK)-XHK-l  ))• 


ISN  0298 
ISN  0299 

ISN  0300 

ISN  0301 
ISN  0302 


ISN  0303 
ISN  0304 
ISN  0305 
ISN  0306 
ISN  0307 
ISN  0308 
ISN  0309 
ISN  03i0 
ISN  0311 
ISN  0312 
ISN  0313 
ISN  0314 
ISN  0315 
ISN  0317 
ISN  0313 
ISN  0319 
ISN  0320 

ISN  0321 
ISN  0322 


ISN  0323 
ISN  0324 
ISN  0325 
ISN  0326 
ISN  0327 
ISN  0328 
ISN  0329 
ISN  0330 
ISN  0331 
ISN  0332 
ISN  0333 
ISN  0334 
ISN  0335 


ISN  0336 
ISN  0337 


ISN  0338 
ISN  0339 
ISN  0340 


• (lJN<LMy,K)  + UN<LMX,K-l  1) 

75  CONTINUE 

USIl.KTEO)  = US< 1 ,KTE0-1 ) + 0.5«( XID-XH KTEO-1 > )• 

• ( UTEU  + UN< 1 .KTEO-1  ) > 

US(LMX,KTEO)  = US{ LMX , KTE 0- 1 ) + 0 . 5*{ X I D-X I < KTEO- 1 ) )* 

• < UTEL  + UN< LMX.KTEO-1  ) ) 

WRITE!  6,203  > 

203  FORMAT! //lOX, 'BLADE  SHAPE',/  8X , ' X I ! K > ' , 8X , ' PCT ' , 6X , ' NU/CA ' , 5X , 

• 'NL/CA'  ,5X, 'H/CA'  ,6X,  'T/CA'  , 6X , ' NU/C  ' , 6X , ' NL/C ' , 6X , 

• 'H/C ,7X, 'T/C' > 

K1  = KLEP  + 1 

OPR  ^ ORT*OMCU 
DO  76  K = K1,KTE0 
US! 1 ,K  > = ORT»US! 1 ,K) 

US!LMX,K)  = ORT-'US!  LMX,K) 

HCA  = !US!1,K)  + US!LMX,K)  )/2. 

TCA  = US! 1 ,K  > - US! LMX,K  ) 

PCT  = XI!K)/X!D 

OU  = US! 1 ,K  )/OPR 

OL  = US!LMX,K>/OPR 

HBC  = HCA/OPR 

TBC  = TCA/OPR 

IF!  :<•  EQ.KTEO  > GO  TO  76 

WRITE!  6,20  4 ) X 1 ! K 1 , PCT , US! 1 , K > , US! LMX , K) , HCA ,TCA ,OU ,OL , HBC ,TBC 

76  CONTINUE 
PCT  = 1 . 

WRITE!  6,204  ) X ID, PCT, US! 1 , KTEO > , US! LMX . KTEO > , HCA, TCA , 

» OU,OL,HBC,TBC 

204  FORMAT! 1PE15.3,0P9F10.4  ) 

RETURN 

C 

C CALCULATE  PERFORMANCE  DATA 
C 

90  CONTINUE 
WRITE!  6, 220  ) 

GMA  = !SB/EM2)  - 1. 

EMF  = 0.5*EM2*!GMA-1  . ) 

EMF  = GMA*EM2/!  TPB*!  1 . +EM2'»EMF  ) > 

DO  91  N = 1,NMX 
W = -DPMI ! N >/! TPB*RO! N ) ) 

U = CONST  - W»RO!N) 

DB  = 57. 296’*ATAN!  ! W-U*RO!  N ) >/!  1 .+RO!  N )*RO(  N > > ) 

SPR  = 1.  + EMF-'DPHnN) 

WRITE! 6,221  > N , ROI N > , DPH 1 ! N > , W , U , DB , SPR 

91  CONTINUE 

220  FORMAT! ///38X,  'PERFORMAi;CE  DATA', 

• ///lOX  , ' N ' ,5X , ' RHO! N ) ' ,6X, ' DPHH N 1 ' ,7X, 'W/U  INF',6X,'U/U  INF', 

• 4X , ' ARCTAN! UN/WO  )',4X,'P02/P01', 

• /68X  , ' ! DEGREES  ) ' ,/  ) 

221  FORMAT! 8X,I3,4X,F7.4,4X,1PE10.3, 

• 5X,0PF8.4,5X,F8,4,6X,F6.2,9X,F6,3> 
RETURN 

C 

C PRINT  MACH  NUMBERS 
C 

96  REL  = SORT!  EM2’*!  1 . +A!  N ) ) > 

CF  = SB/!2.*EM2) 

DO  92  K = 2,KMXM1 
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ISN  0341  DO  93  L = 1 ,LMX 

ISN  0342  US(L.K)  = UB»FAK< K )•(  F ( L , K+ 1 , N > - F(L,K-1,N)) 

ISN  0343  UN(L.K)  = REL*( 1 . +CF*US( L , K > > 

ISN  0344  93  CONTINUE 

ISN  0345  92  CONTINUE 

ISN  0346  WRITE<6.290)  I TR , I TK , N , RO( N > 

ISN  0347  290  FORMAT! 1 H 1 . 3X VALUES  OF  THE  LOCAL  MACH  NUMBER  AFTER  ITR  = ',14, 

• ',  ITK  = ',14,'  AT  RHO(',I2,')  = ’.F7.4) 

ISN  0348  LA  = 1 

ISN  0349  LB  = 20 

ISN  0350  IF(LMX.LT.20)  LB  = LMX 

ISN  0352  97  VIRITE!  6,291  ) ( LOP<  L ) , L = LA,  LB  > 

ISN  0353  291  FORMAT! / 1 X ,' K L=', 12, 1916) 

ISN  0354  WRIT6!6,293) 

ISN  0355  293  FORMAT!/) 

ISN  0356  DO  94  K = 2,KMXM1 

ISN  0357  WRITE!6,Z94)  K , ! UN! L , K ) , L = LA , LB  ) 

ISN  0358  94  CONTINUE 

ISN  0359  294  FORMAT! 1 X , I 2 , 20F 6 . 2 ) 

ISN  0360  LA  = LA  + 20 

ISN  0361  IF! LA.GT.LMX  ) GO  TO  30 

ISN  0363  LB  = LB  + 20 

ISN  0364  IF! LB.GT.LMX  ) LB  = LMX 

ISN  0366  GO  TO  97 

ISN  0367  END 


LEVEL  21.' 


( DEC  72 


) 


OS/360 


FORTRAN  H 


ISN 

ISN 


ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 


ISN 

ISN 

ISN 


COMPILER  OPTIONS 
C 


NAME=  MAIN,OPT=02,L INECNT»60,SIZE=0000K, 

SOURCE, EBCDIC, NOLI  ST, NODECK, LOAD. MAP, NOEDIT, ID, XREF 


0002 

0003 


0004 

0005 

0006 

0007 

0008 

0009 

0010 
0012 

0013 

0014 

0015 


SUBROUTINE  BVALIIBC) 

COMMON  ZT(30),!<I(60),RO(10),F(30,60,10>,A(10>,B(10),D(10>,E(10>. 

* EE(30).FF(30>,SV(30>,DPHI( 10).BV( 10),US(30,60),UN(30,60>, 

• ONDSi  30, 10,2  > ,FAK( 60 ) ,FHV/( 60  > ,UTA( 30  > ,UTB( 30  > ,U7C( 30)  ,BBV<  30), 

» BSA(30),BSB(30),BSC(30  ) , VSA( 30  > , VSBl 30  > , BAV( 30  > ,LOP( 30  ) 

» ,SONIC(10) 

COMMON  DAF .DBF, DCF ,AKLK£ ,RX, OX, CONST 

COMMON  H,DR0,DZT,TDZ.TPB.XIC.XIB,XID,TDR,AKLK,DR02,DZT2 
COMMON  OMCU,RHUB,RTIP.R2.RXE ,RXH,UA,US,TRDZ,ERA,ERB,CKI.KE,SB, 

* BTR.FXB.FXI , EM2 , AMA , AMB , DM , DKA , DKR , DO , DOS , TZ , TTZ , TZS , TRS , TDO 
COMMON  K,N,IDM,  KMX , LMX , NMX , KMXMl , KMXM2 , LMXMl . LMXM2 . NMXMl , 

• NI1XM2,KTE0, IRX,  I TK , ITR , I TRMX . KLEP , I TP R , KUP . KDN 
COMMON  1D(36) 

DO  1 l<  = l.KMX 
IF(XI(K>.GE.XIC)  GO  TO  2 

1 CONTINUE 

2 KLE  = K 

KLEF  = KLE  - 1 
K2  -■=  KTlO  - 1 


C 

C KLEP  IS  THE  LAST  POINT  UPSTREAM  OF  THE  L.  E..  AND  KT£(i  IS  THE  FIRST 
C POINT  DOWNSTREAM  OF  THE  T.  E. 

C 


0016 

0017 

0013 


GO  TO  ; 1 00,200  ) . IBC 
200  DO  20  N = l.NHX 
R2  ’ RO(N)*RO(N) 


C 

C FOR  IBC  = 2,  TAU  = TMAXl R )/C( R > . 

C THE  BLADE  GEOMETRY  USED  HERE  HAS  TMAX  = CONSTANT.  BV(3)  CONTAINS 

C TMAX/CAX,  BV(2)  CONTAINS  SMALL  A,  BV( 1 ) AND  BV< 4 ) CONTAIN  DELTA  CP 

C AT  THE  HUB  AND  TIP,  RESPECTIVELY,  DELTA  CP  ZERO  VARIES  LINEARLY  WITH 

C RHO  BE'iWEEN  THESE  TWO  VALUES. 

C 


ISN 

0019 

TAU  = BV( 3 )/SQRT( 1 . +R2 ) 

ISN 

0020 

AOM  = BV(2)*OMCU 

ISN 

0021 

ACP  = ( eV{  4 )-BV(  1 ) )/l RO( NMX  )-RO(  1 ) ) 

ISN 

0022 

3CP  = ( BV( 1 )-H*BV( 4 > )/( 1 . -H  ) 

ISN 

0023 

DCP  ACP*RO(  N ) + BCP 

ISN 

0024 

D021  K = KLE,K2 

ISN 

0025 

KB  = K - KLEP 

ISN 

0026 

HCP  = DCP/2. 

ISN 

0027 

IF(XI(K).GT.AOM)  GO  TO  22 

1 jN 

0029 

DCP2  = HCP 

ISN 

0030 

DNDS  ( KB  .N, 1 ) = HCP»XI( K ) 

ISN 

0031 

GO  TO  23 

ISN 

0032 

22 

DNDS(KB,N,1  ) = HCP*(AOM+(XI(K)-AOM>»(1.-0.5*(BV(2)+XI<K)/OMCU)) 

• /(  1 .-8V( 2 ) ) > 

ISN 

0033 

DCP2  = HCP*( 1 . -XI< K )/OMCU  )/( 1 .-BV( 2 ) > 

ISN 

0034 

23 

TPS  = 4.*TAU*( 1 .-2.*X1(K)/0MCU) 

ISN 

0035 

0NDS(KB,N,2)  = R2»DCP2/(  1 . +R2  ) + RO( N )*TPS 

ISN 

0036 

DPHKN)  = HCP-*(AOM  + OMCU*(  1 . -BV(  2 ) )/2  . ) 

ISN 

0037 

21 

CONTINUE 

ISN 

0038 

20 

CONTINUE 

ISN 

0039 

WRITE!  6,220  ) BV(  1 ) , BV( 4 > , BV! 2 ) , BV(  3 ) 

ISN 

0040 

220 

FORMAT! //5X, ’BLADE  LOADING  SPECIFIED  FOR  THIS  CASE  HAS  DELTA  CP', 
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* ' =',Fa.4,'  AT  THE  HUB.',F8.4,'  AT  THE  TIP', 

* //lOX, 'LOADING  IS  CONSTANT  OVER  THE  F IRST ' , 2PF6 . 2 , ' PERCENT', 

'*  ' OF  THE  CHORD,  AND  THEN  VARIES  LINEARLY  TO  ZERO  AT  THE  ', 

//lOX,  'TRAILING  EDGE'  , 

//lOX , 'THICKNESS  DISTRIBUTION  IS  THAT  OF  A DOUBLE  PARABOLIC, 

* ' ARC,  V/ITH  TMAX/CAX  = ',0PF8.4,//) 

ISN  0041  GO  TO  10 

ISN  0042  100  CONTINUE 

C 

C THE  BLADE  GEOMETRY  AND  ANGLE  OF  ATTACK  ARE 
C T MAX/CAX  = BV( 1 > + BV(2)/R  + BV(3>*R 
C H MAX/CAX  = BV<4)  + BV<5)/R  + BV<6)*R 
C ALPHA  = BV(7>  + BV<8)/R  + BV(9)*R 
C WHERE  R = RO(N>/RTIP 
C 

ISN  0043  WRITE<6,210>  ( BV{ J J ) , 0 J= 1 , 9 ) 

ISN  0044  210  FORMAT(//5X, 'BLADE  GEOMETRY  AND  ANGLE  OF  ATTACK  SPECIFIED  ARE', 

» /20X,'T  MAX/CAX  = TA  + TB/R  + TC’R,  R - RO(N)/RTIP', 

* /25X, 'WHERE  TA  = ',1PE12.3,'  TB  = ',E12.3,'  TC  --  ',E12.3, 

* //20X,'H  MAX/CAX  = HA  + HB/R  + HC»R', 

/25X, 'WHERE  HA  = ',1PE12.3,'  MB  = ',E12.3,'  HC  = ',E12.3, 

* //20X, 'ALPHA  = AA  + A3/R  + AC*R  (DEGREES)', 

* /25X, 'WHERE  AA  = ',1PE12.3,'  AB  = ',E12.3,'  AC  - ',£12. 3,//) 

ISN  0045  BV(7)  = BV< 7 > / 57 . 29578 

ISN  0046  87(8)  = SV( 8 > /57 . 29578 

ISN  0047  BV<9)  = BV( 9 ) /57 . 29578 

ISN  0048  DO  3 N = 2,NMXM1 

ISN  0049  R = RO(N>/RTIP 

ISN  0050  CBC  = l./SQR7< l.+RO<N)*RO(N)) 

ISN  0051  HMXCA  = BV( 4 ) + 3V(5>/R  + BV(6)*R 

ISN  0052  TMXCA  -■»  BV<  1 ) + BV(2>/R  + BV(3>»R 

ISN  0053  ALPH  = 3V(  7 > + BV(8>/R  + BV(9>'»R 

ISN  0054  DO  4 K = KLE,K2 

ISN  0055  KB  = K - KLEP 

ISN  0056  DNDS<KB,N,1>  = (4. "HMXCA  + 2 . "TMXCA  )*CBC*( 1 . -2. *X I ( K ) /OMCU > - ALPH 

ISN  0057  DNDS<KB,M.2>  = (4. "HMXCA  - 2 . "TMXCA  )"CBC"< 1 . -2 . "X I < K )/OMCU > - ALPH 

ISN  0058  4 CONTINUE 

ISN  0059  3 CONTINUE 

ISN  0060  10  RETURN 

ISN  0061  END 


LEVEL  21.7  ( DEC  72  ) 


OS/360  FORTRAN  H 


ISN 

ISN 


ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 


ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 


COMPILER  OPTIONS 
C 


NAME=  MAIN,OPT=02,LINECNT=60,SIZE=O000K, 

SOURCE .EBCDIC , NOLI  ST, NODECK, LOAD. MAP .NOEDIT, ID. XREF 


0002 

OOOJ 


000-1 

0005 

0006 

0007 

0003 

0009 

0010 
0011 


0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 


SUBROUTINE  RESIDl  ISHO  > 

COMMON  ZT(30),XI(60).RO(10),F(30,60,10),A{10),B<10>,D(10).E(10>, 

* EE(30),FF(30),SV(30).DPHI(10),BV(10),US(30,60),UN(30,60). 

* DNDS(  30, 10,2  > .FAK(  60 ) .FHW( 60 ) .UTA( 30  > .UT3( 30  > .UTC( 30  > .BBV{ 30  > , 

* 3SA<  30  > ,BSB( 30  > ,BSC( 30  > ,VSA( 30 ) , VSB( 30  ) .BAV<  30  ) ,LOP( 30  ) 

* . SONIC!  10  ) 

COMMON  OAF , DBF , DCF .AKLKE , RX , OX .CONST 

COMMON  H , DRO , DZT , TDZ ,TP3,XIC,XIB.XID,T0R , AKLK , DR02 , DZT2 
COMMON  OMCU,RHUB,RTIP.R2.RXE,RXH,UA.UB.TRDZ,ERA,ERB,CKLKE,SB. 

* 3TR ,FXB,FXI,EM2, AMA , AMB , DM , DKA , DKR . DO , DOS , TZ . TTZ , TZS . TRS . TOO 
COMMON  K, N , IDM.  KMX , LMX , NMX . KMXMl , KMXMZ , LMXMl ,LMXM2,NMXM1 , 

* NMXM2,KTE0,  IRX,  ITK  . ITR  , I TRMX  , KL  EP  . I TPP. , KUP  , KDN 
COMMON  ID(36> 

DIMENSION  ERMX( 10),KEP(10),LER<10> 

WR1TE(6.200)  ITR 

200  FORMAT! //5X. 'AFTER  ITR  ^ ',12,'  MAXIMUM  RESIDUALS  ARE', 

* //4X,'N'.7X,' RESIDUAL ',4X,'K'.4X.'L', 

* TX.'AVG  RES',8X.‘AVG  PHI',8X,'AVG  SUM',/) 

DO  1 N = 2,NMXM1 

SUM  = 0. 

SNORM  = 0. 

SBOl  = 0. 

KOUNT  = 0 
ERMX!N>  = 0. 

BTR  = 1.  - EM2*!1.  + RO!N>*RO!N)> 

AKLK  = D!N)*T2S 
DKA  = -B!N)’*TTZ 
ERA  = 1 . + 0 . 5*DRO/RO! N ) 

ERB  = 1.-0 . 6*DRO/RO! M ) 

DO  4 L = 1 ,LMX 
VSA!L)  = 0. 

BSB! L ) = 0. 

BAV!L)  = BTR 
4 CONTINUE 

DO  2 K = 2. KMXMl 
AKLKE  = AKLK/FAKIK) 

DKR  = TRS'-E!  N )/FAK!  K ) 

DO  3 L = 2.LMXM1 
BBV! L > = GAV! L ) 

BAV!L>  = BTR  - SB'*FAK!  K )»!  F ( L . K+ 1 , N )-F  I L , K- 1 , N ) )/T0Q 
VSB! L ) = VSA! L ) 

VSA!L)  = 0. 

IF!BAV!L  l.GT.O.  ) GO  TO  8 
VSA!L)  = 1. 

8 BSC! L > = BSB! L > 

BSB!L>  = FHW!K)«!F!L,K+1,N)-F(L,K,N))  -FHWC K- 1 )•( F ( L . K , N )-F ( L , 

* K-1  .N  ) > 

BSA!L  ) = FHW!K)*F!L,K+1,N>-2.’*FAK!K)*F(L,K.N>  + FHW<K-1 I'FIL.K-I.N) 
TMA  = BAV! L )«!  1 . -VSA!  L ) )*BSA! L ) + BBV! L I'VSB! L >*BSC( L > 

TA  = ABS(TMA) 

TMB  = A! N )•! ! 1 . -VSA! L > )*BSA! L )+VSAl L )*BSB! L > ) 

TB  = ABS!TMB> 

TMC  = -2 .*DKA*! -E! L + 1 ,K-1 .N > + FI L ,K-1 ,N  ) + F! L + 1 ,K,N >-2.»F< L ,K,N> 

* ♦F!L-1,K,N)+F!L.K+1,N)-F!L-1,K+1.N)) 
rc  = ABS! TMC  ) 
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• 


ISN  0047 


ISN 

0048 

TMD  = AKLKE*(F<L+1,K,N)-2.*F(L,K,N>+F(L-1.K,N)) 

ISN 

0049 

TO  = ABS(TMD) 

ISN 

0050 

TME  = DKR*(ERA*F{L,K,N+1  ) + ERB*F(L,K,N-l )-2.*F<L,K,N 

ISN 

0051 

TE  = ABS<  TME  ) 

ISN 

0052 

TOP  = TMA+TMB+TMC+TMD+TME 

ISN 

0053 

BOT  = TA  +TB  +TC  i-TD  +TE 

ISN 

0054 

RES  = ABS(TOP ) 

ISN 

0055 

SUM  = SUM  + RES 

ISN 

0056 

SNORM  = SNOPM  > ABS< F ( L , K , N ) > 

ISN 

0057 

SBOT  = SBOT  + BOT 

ISN 

0058 

KOUNT  = KOUNT  + 1 

ISN 

0059 

IF< ISHO.EQ.O  ) GO  TO  22 

ISN 

0061 

WRITE(6,205)  N.K.L, TMA , TMB , TMC . TMO , TME , TOP . RF S 

ISN 

0062 

205 

F0RMAT<315,1P7E13.3) 

ISN 

0063 

22 

IF<RtS.LT.ERMX(N>)  GO  TO  3 

ISN 

0065 

KER(N>  = K 

ISN 

0066 

LER(N>  = L 

ISN 

0067 

ERMX(N)  = RES 

ISN 

0068 

3 

CONTINUE 

ISM 

0069 

2 

CONTINUE 

ISN 

0070 

SUM  = SUM/FLOAT( KOUNT) 

ISN 

0071 

SNORM  » SNORM/FLOAT( KOUNT) 

ISN 

0072 

SBOT  = SBOT/FLOAT( KOUNT) 

ISN 

0073 

WRITE! 6,210)  N , ERMX ( N ), KER( N ). LERI N ), SUM , SNORM . SBOT 

ISN 

0074 

1 

CONTINUE 

ISN 

0075 

210 

FORMAT! lb,lPE15.3,215,3E15.3) 

ISN 

0076 

10 

RETURN 

ISN 

0077 

END 
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APPENDIX  C 

DICTIONARY  OF  VARIABLES 


ALGEBRAIC 

FORTRAN  SYMBOL 

EQUIVALENT 

DEFINITION,  USE,  COMMENTS 

A(N) 

AE 

See  Equations  (38),  (46) 

ALF 

oc, 

Equation  (16) 

AKL 

K 

Equation  (30) 

AKLK 

/ At 

ip 

AKLKE 

r/ 

AMA,  AMB 

Mach  number  limits  for  which  contours 

1 

are  calculated  if  IOP=4 

B(N) 

-(1  + pV 

BAV(L) 

'K  ' 

BBV 

BE 

See  Equations  (38) , (46) 

BGL 

f’  pA4>  dp 

BKL 

^ B*' 

1 K 

Equation  (30) 

BN 

B 

Number  of  blades 

BSA(L) 

\o-d:)(K  *f"> 

BSB(L) 

BSC(L) 

P^Pk 

BTR 

^ - A 

^90 

BV(I) 

Parameters  for  blade-surface  boundary 

conditions 

CALT 

Co./ i-r 

CE 

See  Equations  (38),  (46) 
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FORTR.AN  SYMBOL 

ALOE BRA  1C 
EQUIVALENT 

DEFINITION,  USE,  COMMENTS 

CKL 

C*" 

'-K 

Equation  (30) 

CKLKE 

= RKLKE 

CONST 

C 

See  Equation  (10) 

D(N) 

(1  V 

DB 

Used  in  Subroutine  OUTPUT 

DAF.[)BF,DCF 

See  Equation  (52) 

DKA 

(i+p^)  at/iat; 

DKL 

dI; 

Equation  (30) 

DKR 

J +p'‘  / At 

DM 

i 

j 

Mach  number  interval  for  which  contours 

are  calculated  if  IOP=4 

DNDS(K,N,J) 

Surface  slopes  for  IBC=1,  loading  and 

thickness  parameters  for  IBC=2,  see 

Equations  (17)  and  (18) 

DPMI (N) 

A(P  (p) 

DQ 

At 

1 

DQS 

ti 

<1 

DRO 

Ap 

DR02 

(Apf 

DZT 

Al^ 

DZT2 

E(N) 

1 + 

E(N)  is  defined  to  be  zero  for  IDM=2 

FIMTC. 

EMX 

EM2 

M ^ 

' OO 

ERA 

1+  Ap/zp 

C-2 


FORTRAN  SYMBOL 

ALGEBRAIC 

EQUIVALENT 

DEFINITION,  USE,  COMMENTS 

ERB 

/ - ^p/ zp 

F(L,K,N) 

4>CK,z,p) 

FAK(K) 

FHW(K) 

^ K 

tJ.X  »>.X  l>J,X 

FPA,FPB,FPC 

FXB 

Used  to  locate  upstream  and  downstream 

FXI 

edges  of  the  grid  - see  Equation  (16) 

GMA 

n 

H 

^ - PJPr 

IBC 

Indicator  for  blade-surface  boundary 
conditions 

I BN 

8 

I DM 

=2,  3,  for  two,  three-dimensional 
calculations 

lOP 

Used  to  select  output  option 

IRPT 

Print  control 

IRX 

Equals  2 or  1 if  there  are  or  are  not 
any  hyperbolic  points  on  the  line  being 
relaxed 

IRXP 

1 

Controls  intervals  at  which  residuals 
are  printed 

IRXPT 

Print  control 

ISAVE 

= 1,  0 if  results  are,  are  not  to  be  saved 
on  tape 

ISHO 

= 1 if  residuals  are  to  be  shown  at 
every  grid  point 

I START 

Selects  start  option 

ITK 

Counter  for  iteration  in  Z -direction 

ITKMX 

Maximum  number  of  iterations  inz- 
direction 

( 

:-3 

■ORTRAN  SYMBOL 


ITRMX 


ALCLBRAIC 
LQU I VALl-NT 


I)I-.FINITION,  list;,  COMMENTS 

= 1 if  elliptic  relaxation  factor  is 
tapered,  see  Equation  (A-4) 

('ounter  for  iterations  in  p -direction 

Maximum  number  of  iterations  in  p -direc- 
t ion 


Print  control 


KDN.KUP 


! Downstream,  upstream  limits  within  which 
Mach  number  contours  are  calculated  if 
10P=4 


KMXMl 


KMX  - J 


Number  of  grid  points  in  the  Z -direction 


KMXM2 


KMX  '2 


KW,KWMX 


Last  K-station  upstream  of  leading  edge 

First  K-station  downstream  of  trailing 
edge 

Counters  for  locations  of  Mach  number 
contours 


LMXMl 


i-MX  - 1 


Number  of  grid  points  in  the  ^ -direction 


LMXM2 


1-MX  -Z 


NMXMl 


NMXM2 


/VMX-  Z. 


Number  of  grid  points  in  the  ^-direction 


/ - OBV 


Print  control 


Print  control 


Relaxation  factor  for  circulation,  see 
Equation  (53) 
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ALGEBRAIC 

FORTRAN  SYMBOL 

EQUIVALENT 

DEFINITION,  USE,  COMMENTS 

OMCU 

/ iX^ 

Value  of  Z at  trailing  edge  of  blades 

OPR 

OPT 

See  Equation  (A-1) 

ORT 

OX 

/ - }/Ca) 

Q 

Z 

REL 

Inlet  relative  Mach  number  at  r" 

RELSQ 

(.  /cu„f 

RES 

"/?K 

See  Equation  (73) 

RHUB 

RO(N) 

/O  Ljr/ix„ 

RTIP 

Pt  = 

RX 

iA> 

Relaxation  factor,  see  Equation  (33) 

RXE 

Relaxation  factors  used  at  elliptic. 

RXH 

hyperbolic  points 

RXM 

1 -60 

R2 

SB 

(2(4  f;  mI 

SPR 

/ i’o, 

SV(L) 

Temporary  storage  for  (j) 

TA-TE 

See  Equations  (68) -(72) 

TALF 

X Oc, 

See  Equation  (16) 

TDQ 

ZAt 

TDR 

ZA  p 

TDZ 

C-5 


FORTRAN  SYMBOL 

ALGEBRAIC 

EQUIVALENT 

DEFINITION,  USE,  COMMENTS 

TMA-TME 

Absolute  values  of  TA  through  TE 

TPB 

xrx  / B = 2;, 

TP  I 

zn 

TRS 

(At  / Ap)^ 

TRDZ 

XpA  TT 

TTZ 

At  / 2 A/: 

TZ 

At /At; 

' 

TZS 

CAt/At;f 

U 

Used  in  Subroutine  OUTPUT 

UA 

P 

UB 

[2AT  O + p’;] 

ULEL 

1 

ULEU 

i 

UN 

\ 

UR 

Tr/  tX„ 

US 

LXg/  Wo 

i 

USCRIT 

Sonic  value  of 

UTA.UTB.UTC 

i 

f/.i.  hJ  ,L.  ^I.L. 

Values  of  <?«  , t«.i  > from 

UTEL 

UTFU 

previous  Z -iteration  • 

t 

VSA 

M-k 

VSB 

W 

Used  in  Subroutine  OUTPUT 

XI  (K) 

2 = U)X/tX„ 

XIB 

Zb 

See  Equation  (16) 
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FORTRAN  SYMBOL 


ALGEBRAIC 

EQUIVALENT 


DEFINITION,  USE,  COMMENTS 


C-7 


